三维“体点”散热问题建模仿真分析

徐仲强，王远弟

(上海大学数学系，上海 200444)

1 引言

随着电子元器件尺寸越来越小，功率也不断攀升，使得元器件散热变得愈发重要，元器件的传热优化问题受到广泛关注[1-3]。

现在这个问题的有效解决办法之一是在电子元器件内部构造高传热材料传导路径，通过插入高传热性材料(如金刚石或碳纤维)，可以更有效地排出内部的热量。如何利用有限数量的高传热性材料来构造热传导路径，将平均温度最小化转化为一个优化问题，在热传导优化中，称为“体点”问题(VP)。“体点”问题是热传导优化的基础问题，最初由Bejan定义。

Bejan等人于1996年[4]提出了基于构型理论的树状分叉网络，讨论了其在微电子器件冷却散热优化方面的应用，对平面体点问题进行了优化。他的基本思路是从最小的二维面积单元开始(通常与工艺上能达到的最小制造单元相当)，对该面积单元利用优化方法使其热阻最小，达到其最优外形；然后利用优化方法对这个最小的面积单元进行第一次组合，得到第一次组合体结果，从第二次组合体开始，新的组合体一定包含两个较小的两个优化后的组合体。按照这样的优化方法进行下去，直到经过若干次的组合优化后得到的最优组合体能覆盖住所给的面积。Bejan的构形理论建立在主次干道相互垂直的基础上，并对每一干道的宽度和长度进行优化，得到了在内热源均匀、高传热材料和基体材料的传热系数比值较大的情况下高导材料最优布置方案。

过增元等人以提高传热效率为优化目标，基于生命演化的自然原理提出并发展了仿生优化方法[7]，定义了火积(描述热传导能力的物理量)[5]、热量传递势容和热量传递势容耗散函数[6]，提出了最小热量传递势容耗散原理，得出了当传热系数最佳分布时，全场的温度梯度应处处相等。仿生优化方法所遵循的梯度均匀化原则的理论推导利用了传热系数的连续性[6]，对内热源均匀且高传热材料和基体材料的传热系数比值较小、非均匀内热源[8]、具有相变的非稳态状态[9]等问题都有较好的优化结果。

2018年，王远弟和张俊顶等人[10]在讨论二维传热问题时，以全场温度均值最低为优化目标，提出并证明了高传热系数材料最优分布时的曲面面积极小化原则，最低全场温度均值时高传热材料布置对应最小温度场曲面面积，利用数值模拟说明了极小曲面法的合理有效性，并对比了极小曲面法和仿生优化方法在不同条件下的优化效果。

受二维仿生优化以及极小曲面启发，本文的优化方向分别是提升传热效率和最小化全场温度均值，对于三维“体点”问题，得出了三维高导材料填充标准。在此基础上对两种方法填充，对比剖析两种方法模拟出的构造，可知三维模型的填充结构与二维类似，受到高传热材料传热比以及填充量的影响，得到的三维模型高传热材料结构对于实际问题具有借鉴意义。

2 数学模型

底部中心模型如图1所示。

图1 底部中心模型

图1左图模型a和右图模型b描述了一个空间元器件材料“体点”传热问题。立方体边长为L，内均匀填充热源q。高导材料和初始材料的传热系数分别为Kp和K0，高导材料填充量为Vp，填充率为η(Vp/V)。模型a底部有长度为δ(δ≪L)的正方形开口，边界绝热，开口温度恒定U0。模型b中心镂空，镂空处保持恒定温度U0，边界绝热。 求区域内高传热系数材料最优填充，使得全场温度均值最低。一直处在运行中的器件，可以视为稳态，所以考虑平衡态下的传热问题。设u=u(x，y，z)为温度函数，为了最小化全场温度均值。记

(1)

为全场温度均值。其中|Ω|为空间区域Ω的度量即体积，易知|Ω|恒定，省略它对结果无影响。于是三维“体点”传热问题的数学模型a、b可化简为

(2)

在边界面各个方向上，对于底部开口的模型a：

(3)

对于中间镂空的模型b

(4)

其中Γ11、Γ21、Γ22、Γ23分别为模型a的开口、前后绝缘壁、左右绝缘壁和上下绝缘壁，Γ12、Γ24、Γ25、Γ26分别为模型b的开口、前后绝缘壁、左右绝缘壁和上下绝缘壁，式 (2)为优化目标，式(3)(4)分别为模型a、b的边界条件。取传热口处不变温度为u=0，否则做平移变换即可得。

3 理论分析

这里讨论两种三维模型填充方法，即仿生优化和极小曲面方法。注意到模型a和b仅仅只是在中间温度的区别，外部边界条件类似，所以这里以模型a为例作理论推导。

3.1 极小曲面方法

数学上，把R3中平均曲率为零的曲面称为极小曲面，该函数u=u(x，y)满足偏微分方程

(5)

相应地，三维空间上的函数u=u(x，y)满足极小曲面方程

(6)

(7)

(8)

计算导数

(9)

其中v={cosα，cosβ，cosγ}是∂Ω的外法向量，由w∈W的任意性，可得

(10)

由式(10)不难计算出当全场温度均值最低时，温度函数满足极小曲面方程。

(11)

3.2 仿生优化方法

对于三维VP问题，同样引入了火积[6]

Zdis=∭k|∇u|2dA

(12)

其中，Zdis表示热量传递过程中的火积耗散，k为传热系数，u=u(x，y，z)是传热区域Ω上平衡态下的温度函数， ∇u为温度梯度。火积耗散越小，则温度梯度场越均匀，给定区域的温度也越低，因而传热优化程度就越高。因此，在对元器件模型作网块划分后，高传热材料应该首先放在火积最大的微元位置，这些点相当于火积最“突出”的点，当把高导材料填充在这些位置点时，必然最大程度地减少整体火积。

3.2.1 三维球状元件的散热问题

最小热量传递势容耗散原理指出传热系数或密度为最佳径向分布时，热量传递势容耗散最小，即满足

(14)

3.2.2 三维体点散热仿生优化

设u=u(x，y，z)是传热区域Ω上的温度函数，以火积耗散最小为优化目标，只需分析min∭Ωk|∇u|2dV，这里的|Ω|指空间区域Ω 的度量即体积。满足第2节中模型a的边界条件。这里dV是体积微元，u=u(x，y，z)，u=u(x，y，z)∈

kuxx+kuyy+kuzz+kxux+kyuy+kzuz=0

(15)

当k为常数时，温度函数u=u(x，y，z)满足拉普拉斯方程。

4 高传热系数材料优化准则

4.1 极小曲面

根据3.1的理论推导可知，最低全场温度均值时温度函数图像对应的超曲面是“极小曲面”。因而，三维超曲面极小化原则指明高导材料最优布置应该使得温度函数图像的“面积”趋于极小化。由于每一块高导材料的放置都会影响先前的温度场，温度函数图像更新。故极小曲面原则即要求不断地极小化温度函数图像的“面积”。

然而温度函数图像“面积”极小化难以量化，高导材料的最优布置只通过极小曲面原则是无法明确的，所以需要进一步简化超曲面“面积”极小化原则。

数值模拟过程中利用含有源项的稳态热传导方程

∇(K(x，y，z)∇u)+q=0

(16)

设内热源为q=10000W/m3，图2左图是初始温度场，其中邻接网块点间温差小于1K。图2右图是填充高导材料后的温度场，存在温度为338.5959K的网块点，其邻接上、下、左、右、前、后六个网格点上的温度分别是338.9062K，338.3786K，338.6786K，338.4355K，338.6786K和338.4355K。即ux<1，uy<1，uz<1所以有uxuyuxy≪1，uxuzuxz≪1，uyuzuyz≪1。这是邻接网块点间温差小于的一个例子，填充后温度场各网块均符合此现象。易知填充高导材料后温差变小，式(11)表达的极小曲面方程可近似看作调和方程

uxx+uyy+uzz=0

(17)

图2 填充前模型a温度场(左图)

模型a极小曲面法传热比为 300，填充率为8%温度场(右图)

通过将最低均值温度时的温度函数近似为调和函数，便可利用调和函数的平均值定理[11]，即

(18)

上式为调和函数平均值定理满足的必要条件，记

(19)

4.2 仿生优化

仿生优化方法[7]将高导材料的布置分为进化与退化。在温度梯度最大位置进化，梯度最小处退化淘汰，温度梯度趋向均匀化。对(12)式作离散化处理，对立方体元器件作步长为n的分割后，目标函数转化为

(20)

图3 (i，j，z)微元方块与其右侧方块界面

对于微元方块(i，j，z)处，其梯度模的二次方

(21)

对于内部微元方块，比如图3中微元方块(i，j，z)，考察其右界面，根据傅里叶定律(Q=-kAdT/dx)有

(22)

即

(23)

同理其它非边缘界面上相应温度

对于边缘界面，比如右边缘微元方块，其右边界温度Tw=T(i，j，z)。

5 数值计算

5.1 数值计算细化

因为模型a和b区别仅仅在于开口构造不同，数值计算本质上原理趋同，不失一般性，这里以模型a为例阐明数值模拟细则。取基体材料传热系数处处为1，开口边长占元器件边长的1/5，开口温度为300K。

极小曲面方法根据式(19)，采用中心点周围相邻一层的球体邻域上的积分平均值来进行数值模拟。数值模拟 (参考图19) 可知，按球体邻域积分均值所得全场温度均值效果优于正方体邻域，这与球体邻域相邻方块更能代表中心体温度信息有关，正方体邻域中，四角相对较远的方块会冲淡相邻方块含有的中心体温度变化信息，从而使得正方体邻域效果相对不如球体邻域，所以后面以球体邻域进行计算。

填充过程中，利用对称性，空间上每个高导材料填充点有四个对称点。于是，三维“体点”传热问题的数值模拟过程可细化为：

1)预定义参数，包括传热系数，网块步长；

2)根据(16)式求温度场散布；

3)利用(18)式求得中心微元体与其球体邻域温度均值的差值(绝对值)场；

4)在差值最大处放置高导材料，更新传热系数值；

5)重新利用步骤2)更新温度场，当达到填充比例η时终止填充。

仿生优化方法中，整个模拟计算过程具体操作步骤不同之处在于上述步骤3) 4)，仿生优化方法中，特殊的操作步骤为：

3′)根据式(20)计算每个微元体的火积，得到火积场；

4′)根据放置原则，由火积场选择火积最大的微元块作为高传热系数材料最优填充位置，并改变选取位置的传热系数值。

5.2 优化结果

本文是基于模型a和b采用极小曲面方法和仿生优化方法的数值模拟。三维模型相比二维模型，计算量激增，具体到式(11)求解温度场的线性方程组系数矩阵从二维的n2×n2跃升为三维的n3×n3维度(其中n表示每一边长剖分后的网格/块数)。

下文分析中，导热系数比值较大以高导热系数材料与基体材料导热系数比值等于300为例，导热系数比值较小以导热系数比值等于4为例，填充率取8%。

5.2.1 模型a

模型a极小曲面

传热系数比值比较大时，图2左图表示初始温度场，温度最高值为344.6795K，在上方四角，全场温度均值为T=340.4682K。图2右图表示极小曲面原则优化后温度场，全场温度均值为T=301.0257K，同比下降39.4425K；最高温度为301.5658K，同比下降43.1137K。可见高导材料放置后，最高温度和全场温度均值都有明显降低。

图4 模型a极小曲面法传热比为 300，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

传热系数比值比较小时：优化结果如图6所示，图5为相应的温度场散布图像。优化后全场温度均值T=314.0793K，传热区域内最高温度为317.654K。

模型a仿生优化

传热系数比值比较大时：梯度均匀化原则优化结果如图7所示，优化后全场温度均值仅为T=301.0739K，传热区域内最高温度为301.7303K。

图5 模型a极小曲面法传热比为 4，填充率为 8%温度场

图6 模型a极小曲面法传热比为 4，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

图7 模型a仿生优化法传热比为 300，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

传热系数比值比较小时：优化结果如图8所示，优化后全场温度均值T=313.9304K，最高温度为317.5472K。

图8 模型a仿生优化法传热比为 4，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

5.2.2 模型b

模型b极小曲面

传热系数比值比较大时 图9左图是模型b填充前的温度场，由于模型b中间镂空开口温度恒定，故未填充之前Z 轴方向上每一层XY平面温度相等，且从内向外温度越来越高，模型b最高温度出现在到散热竖直开口距离最远的四角处，大小为305.6639K，全场温度均值为T=304.3392K，注意到模型a优化前温度最高处大小为344.6795K，全场温度均值为T=340.4682K，因为模型b镂空开口比较大，所以比模型a的散热效果要好很多，在不填充传热材料时也能达到较低的体平均温度；另外，由于模型b未填充前每一层温度完全相同，这时极小曲面方法会有短暂的失效，这时根据其它填充经验，第一个点填充在开口附近便能得到比较理想的填充效果，于是便采用第一点(准确地说是根据对称性预填充四个点)预填充人为得改变初始温度场，这时候极小曲面方法便可在此基础上达到最终填充完整效果。填充效果如图10所示，图9右图是模型b填充后的温度场。根据极小曲面原则优化后全场温度均值仅为T=301.5272K，传热区域内最高温度为301.9917K。

传热系数比值比较小时：优化结果如图11所示，图12是模型b填充后的温度场。根据极小曲面优原则优化后全场温度均值为T=303.3659K，最高温度为304.6382K。通过填充效果图发现，这时候的填充效果并不是十分理想，在后面的分析中，会对比仿生优化方法说明模型b极小曲面方法局限性。

图9 模型b填充前温度场(左图)、模型b极小曲面法传热比为 300，填充率为 8%温度场(右图)

图10 模型b极小曲面法传热比为 300，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

图11 模型b极小曲面法传热比为 4，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

图12 模型b极小曲面法传热比为 4，填充率为 8%温度场

模型b仿生优化

传热系数比值比较大时 优化结果如图13所示，优化后全场温度均值T=301.1784K，最高温度为301.6163K。

图13 模型b仿生优化法传热比为 300，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

图14 模型b仿生优化法传热比为 4，填充率为 8%三视图：主视图(左上)、左视图(右上)、俯视图(左下)、三维图(右下)

传热系数比值比较小时 优化结果如图14所示，优化后全场温度均值T=302.8593K，最高温度为304.0459K。

5.3 结果分析

以上仅仅对部分结果进行了展示，本节重点对结果进行分析：

5.3.1 两种方法对比分析

在5.2节中，展示了两个模型两种方法的填充率为8%，传热比分别为300 和4 时候的情况，对于其它传热系数比值以及填充比下两个模型两种方法的模拟对比见下文分析。

优化效果(体平均温度变化曲线)

模型a：图15～图18分别描述了传热比在2～400，1～126，127～269，270～400，填充率在3%～8%时两种方法的温度均值差，这里横轴为传热比，纵轴为温度均值差。虚线以下部分表示极小曲面方法所得温度均值更低。

图15 传热比1～400，填充率3%～8%下两种方法计算的温度均值差值曲线

图16 传热比1～126，填充率3%～8%下两种方法计算的温度均值差值曲线

图17 传热比127～269，填充率3%～8%下两种方法计算的温度均值差值曲线

图18 传热比270～400，填充率3%～8%下两种方法计算的温度均值差值曲线

从图15～图18四个图中的计算结果可以看出，不同填充条件下，极小曲面和仿生优化方法互有所长：①图15表明优化曲线稳定后，传热比值在170～400，填充率为5%～8%时，极小曲面原则效果较好。②如图16所示，当传热比值偏小，大约低于126K时仿生优化效果更佳。③图17和图18表明，除了当传热系数比值较大(127K以上)时填充率为3%和4%时整体上仿生优化效果好；填充率为5%～8% 时，传热比127～269时，两种方法优化效果基本相同，传热比270～400时，极小曲面方法优化效果优于仿生优化方法，两者的最大全场温度均值差值在0.1K以内。

模型b：图19表示传热系数比值在2～400之间，填充率为3%～8%时极小曲面方法采用球形邻域与方形邻域的温度均值差值曲线，显然，传热系数比值在50以上时，球形邻域优化效果要比方形邻域好，最大差值达到了0.5K，所以在取平均值时选择的邻域为球形邻域；图20表示传热系数比值在2～400之间，填充率为3%～8%时极小曲面方法采用球形邻域极小曲面方法与仿生优化方法的温度均值差值曲线，由图可见，模型b中，仿生优化效果优于极小曲面方法，且随着高导材料比值越大，两者差异越小。

图19 球形邻域方形邻域温度均值差值曲线

图20 球形邻域仿生优化温度均值差值曲线

公共填充率

在本节中，把两种方法优化结果共同填充点在总填充点中的比例作为公共填充率，公共填充率描述了两种方法共同填充点，在某种程度上，可以作为填充优化的置信度。

模型a：图21表示传热系数比值在1～400之间，填充率为3%～8%时极小曲面方法和仿生优化方法优化结果中填充共同点的比例，由图可见，当传热系数大于一定程度时，共同填充率趋于平稳；图22表示传热系数比值在300～400之间，填充率为3%～8%时极小曲面方法和仿生优化方法优化结果中填充共同点的比例，共同填充率在30%～65%之间，且填充率为5%时共同填充率相对最高。

图21 模型a传热比1～400，填充率3%～8%两种方法填充共同点比例

图22 模型a传热比300～400，填充率3%～8%两种方法填充共同点比例

模型b 图23表示传热系数比值在1～400，填充率为3%～8%时极小曲面方法和仿生优化方法优化结果中填充共同点比例，由图可见，整体上共同填充率比较平稳，填充率为3%时，两种方法会间断出现无共同填充点现象，这与填充比较低时填充点少有一定关系；图24表示传热系数比值在300～400，填充率为3%～8%时极小曲面方法和仿生优化方法优化结果中填充共同点比例在3%～40%之间，且填充率为7%时共同填充率相对最高。

图23 模型b传热比1～400，填充率3%～8%两种方法填充共同点比例

图24 模型b传热比300～400，填充率3%～8%两种方法填充共同点比例

由公共填充率分析可知两个模型两种方法优化结果的共同填充率比较平稳，整体上在20%以上，结合两种方法的优化效果，可以相信两种方法都能达到比较理想的优化效果。

5.3.2 20-50维度对比分析

图25左右图分别表示模型a传热系数比为18，填充率为3%时极小曲面方法和仿生优化方法方法的填充效果，此时，两种填充微元均集中在开口附近，这与20 维度三维效果一致。

图25 模型a传热比为 18，填充率为3%极小曲面法填充(左图)、仿生优化法填充(右图)

图26左图表示模型a传热系数比值为300，填充率为5%时极小曲面方法的填充效果，图26右图表示模型b传热系数比值为150，填充率为5%时极小曲面方法的填充效果，此时，两种填充微元均比较分散，呈细长分布状，这与20 维度三维效果一致。

图26 模型a传热比为300，填充率为5%极小曲面法填充(左图)、模型b传热比为150，填充率为5%极小曲面法填充(右图)

5.3.3 2D-3D对比分析

图27 二维传热系数比为300，填充率为8%极小曲面法填充(左图)、仿生优化法填充(右图)

图27左图和右图分别表示二维传热系数比值为300，填充率为8%时极小曲面和仿生优化方法填充，对比图4以及图7发现，传热比值偏大情况下，三维空间上高导材料有类似二维优化结果的长条，这与高传热系数材料能显著改变区域内的传热效果有关，长条形态使热量沿着高传热长条更快得流向散热口处。

图28正上图、左下图、右下图分别表示三维模型a传热比为300，填充率为8% 极小曲面法50 维度填充结果的Z=1半横切面、X=23纵切面以及Y=3纵切面，可见，在三维空间元器件的二维切面上，也有类似二维优化结果的细长分布长条，这在数值上说明三维极小曲面和仿生优化理论与二维结论有一定相似性。

图28 模型a极小曲面法传热比为300，填充率为8%，Z=1半横切面(正上图)、X=23纵切面(左下图)、Y=3纵切面(右下图)

图29 二维传热系数比值为4，填充率为8%极小曲面法填充(左图)、仿生优化法填充(右图)

图29左图和右图分别表示二维传热系数比值为4，填充率为8%时极小曲面和仿生优化方法填充，对比图6以及图8发现，当传热系数比值较小时，三维空间上高传热系数材料有类似二维优化效果，高传热材料集中在开口附近。

图30 模型a极小曲面法传热比为4，填充率为8%，Z=2半横切面(正上图)、X=18纵切面(左下图)、Y=25纵切面(右下图)

图30正上图、左下图、右下图分别表示三维模型a传热比为4，填充率为8%极小曲面法50 维度填充结果的Z=2半横切面、X=18纵切面以及Y=25纵切面，可见，在三维空间元器件的二维切面上，也有类似二维优化结果的开口集中形状，直观得说明了数值上三维极小曲面和仿生优化理论与二维结果有一定相似性。

6 结论

本文借鉴二维“体点”模型研究方法，对三维“体点”问题中高导材料的最佳填充问题进行研究，根据三维极小曲面方法和三维仿生优化理论找到了三维填充优化准则。数值计算结果表明，三维“体点”传热问题中高导材料填充问题可以通过极小曲面和仿生优化方法加以解决，有如下结论：

1)三维“体点”最低全场温度均值时高导材料分布对应最小温度场曲面面积，即极小曲面的面积。

2)三维“体点”传热问题中高导材料的填充准则与二维最优分布填充准则一致，可根据超曲面面积极小化原则(极小曲面原则)以及火积耗散最小原则(仿生优化原则)来确定。

3)三维“体点”传热问题中高传热系数材料最优区域分布与二维“体点”传热问题中高传热系数材料最优区域分布具有很大的相似性。从三维立体结构规律上看，极小曲面和仿生优化方法得到的高导材料最优填充相像。传热比值偏大时高导材料趋向长条状；传热系数比值偏小时高导材料开口处聚集。从三维“体点”的二维切面上看，高传热系数材料最优区域分布中存在一些切面与二维“体点”传热问题中高传热系数材料最优区域分布具有很大的相似性，满足二维分布规律。

4)通过对比可知，三维“体点”传热问题中超曲面面积极小化方法和仿生优化方法对“体点”传热问题的最优分布有一定比例的共同点，且两种方法的优化效果各有优势。当开口竖直镂空时，仿生优化方法优化效果好；当开口仅存在底部时，填充率为5%～8%，高传热系数比值在127～400时极小曲面方法优化效果较好。