微纳尺度体点导热的拓扑优化
李含灵; 曹炳阳
正则化几何因子
【期刊名称】《《物理学报》》
【年(卷),期】2019(068)020
【总页数】11页(P11-21)
【关键词】体点导热; 拓扑优化; 微纳尺度; 声子玻尔兹曼输运方程
【作 者】李含灵; 曹炳阳
【作者单位】清华大学工程力学系  热科学与动力工程教育部重点实验室  北京 100084
【正文语种】中 文
1 引 言
以芯片为代表的电子器件,在当代社会中扮演着不可替代的重要角.随着电子器件向小型化、集成化的方向发展,器件的特征尺寸已经降低至微纳米量级[1],器件内的功率密度大幅增加.2012年,芯片内的平均热流密度便已达到250 W/cm2左右[2].高功率密度会让电子器件温度升高,产生局部热点,导致其可靠性和寿命显著下降[3].因此,电子器件的热管理已经成为相关领域发展面临的关键挑战之一[4-7].根据电子器件内热点分布和传热结构布局等因素,进行散热仿真和优化改进,对实际产品的设计和制造有着重要意义.
在芯片内部填充高导热材料以加强导热是提高芯片散热能力的主要途径,考虑到芯片内空间宝贵,高导热材料的添加数量必须受到限制,这样的芯片散热问题可以被抽象为体点(volume-to-point,VP)导热问题[8],如图1所示.给定体积的区域代表芯片,区域内的均匀内热源代表焦耳加热,产生的热量通过边界处一小段温度为 T0 的低温热沉流出,其余边界绝热.现提供数量有限的高导热填充材料,希望寻合适的填充材料分布方式,构造高导热通道,让整个系统的温度(平均温度或最高温度)降到最低.关于体点导热问题,人们采用多种优化方法进行了研究,包括构型理论[8]、仿生优化[9]、拓扑优化[10-13]、模拟退火和遗传算法[14]等.尽管这些方法都能有效降低系统温度、提高散热能力,但所得的高导热填充材料分布方式不尽相同,这与优化问题自身带有的局部最优性和多解性有关.当优化问题改变时,不同方法的性能优劣也可能发生变化.
在各类优化方法中,拓扑优化方法因其能够改变结构的拓扑构型、具备理论上的最高设计自由度而受到人们的青睐.利用拓扑优化,设计者不仅可以实现设计对象某方面性能的显著提升,而且所花费的设计时间更少,能得到传统经验之外的具有启发性的新结果[15].体点导热问题的拓扑优化会得到填充材料充分伸入内部区域、整体呈树状的材料分布,其散热性能明显优于其他方法所预测的结构[11].
图1 体点导热问题示意图Fig.1.The schematic diagram of the VP problem.
已有的导热拓扑优化方法,都建立在经典的傅里叶导热定律的基础上[13],这在宏观尺度下是合理的,但对微纳米尺度下的导热过程,傅里叶导热定律不再成立[16-22],热仿真及热优化必须考虑非傅里叶效应.声子是电子器件半导体中的主要载热子[23],声子弹道输运和声子相干是微纳米尺度下发生非傅里叶导热的两个重要原因[5].声子弹道输运的强度可由努森数 Kn 来表征,其定义为材料声子平均自由程l和系统特征尺寸L的比值,即Kn=l/L.宏观尺度下,系统尺寸远大于声子平均自由程(Kn→0),声子间散射充分,声子按扩散方式传递,导热过程遵循傅里叶导热定律.但当系统尺寸与声子平均自由程相当(Kn～1)甚至更小(Kn＜1)时,声子不经历散射而直接抵达边界,这样的过程称为弹道输运,会导致傅里叶导热定律失效[24].在微纳米结构中,部分声
子发生弹道输运,其余声子仍按扩散方式运动,对应的导热机制称为弹道-扩散导热[25,26].声子弹道输运会引起热导率的尺寸效应,表现为系统的等效热导率与系统尺寸正相关,尺寸减小会造成等效热导率降低[27,28],系统热点温度升高[29].此外,在给定温度边界处,由于弹道输运的非平衡本质,边界处的温度并不连续,会出现温度跳跃现象[30,31].声子相干则主要出现在特征尺寸和声子波长相当的系统中[32],例如石墨烯这样的二维材料和超晶格这样的特殊结构[33-35].声子相干通过改变声子速度、弛豫时间和态密度等因素来影响导热过程[36].对常用的半导体材料硅,室温下声子平均自由程约为300 nm[37],声子波长小于10 nm[17],而目前大多数电子器件的特征尺寸在100 nm量级[1],故本文主要关注声子弹道输运引起的非傅里叶导热.关于微纳米尺度下弹道-扩散导热过程的优化,前人曾结合动力学理论和拓扑优化来处理[38],但其优化目标只是特定点之间的温差,无法考虑系统整体的温度.目前,尚缺少适用于微纳米尺度下弹道-扩散导热过程的一般性优化方法,关于弹道效应对优化结果的影响也缺乏系统的认识.
本文结合声子玻尔兹曼输运方程(Boltzmann transport equation,BTE)和固体各向同性材料惩罚(solid isotropic material with penalization,SIMP)方法,发展了微纳尺度下弹道-扩散导热的拓扑优化方法.用声子BTE代替傅里叶导热定律作为导热本构方程,以反映声子弹道输运; 用SIMP方法对声子平均自由程的倒数进行插值以引入拓扑优化,并通过对相对密度变量全局梯
度施加显示约束的方法来确保解的适定性和网格无关性.将此方法用于体点导热优化,发现弹道-扩散导热机制下,拓扑优化预测的填充材料分布明显不同于傅里叶定律下经典的树状分布,且会随 Kn 的变化而变化,这与SIMP方法的插值方式和声子弹道输运有关.
2 优化方法
2.1 声子玻尔兹曼输运方程
声子的运动规律符合声子BTE,稳态、弛豫时间近似下的声子BTE写作
其中f是声子分布函数; vg 是速度矢量; f0 是平衡态声子分布函数,即玻-爱因斯坦分布; τ 是弛豫时间.联合求解(1)式和能量守恒方程,就可以获得声子的微观运动规律和分布函数,进一步可计算得到温度、热导率等宏观热性质.
本文选择离散坐标法(discrete ordinate method,DOM)[39-41]来求解声子BTE,此方法最早用于辐射输运方程的求解,相关计算方法的发展较为成熟.根据声子和光子的类比,可定义声子辐射强度其中ω表示声子频率,D(ω)表示态密度,“ ”表示对声子极化的求和.利用此定义,可以将声子BTE转化为声子的辐射输运方程(equation of phonon radiative transfer,EPRT)[42]
其中s是单位方向矢量,为声子消光系数,是声子平衡态辐射强度.DOM解EPRT的基本思路是[18]：假设I在某个立体角单元内不随方向变化,则整个立体角空间可被划分为多个控制角度,每个控制角度对应的单位方向矢量是 si ,每个离散方向上的EPRT可写作si·∇Ii=κ(I0-Ii),求解此式∫可得不同方向的辐射强度 Ii ,再利用加权求和得到其中 dΩ 表示立体角的微分,wi 是与角度划分方式相对应的权系数,N是控制角度的总数.稳态下,含内热源的能量守恒方程为 其中 q 是热流密度矢量,是单位体积的内热源功率.利用热流密度和声子辐射强度的关系,最终得到DOM求解EPRT的导热方程组为
(3)式中包含N+1个等式,与未知的声子辐射强度数量相对应,方程组封闭.
得到声子辐射强度后便可计算温度.声子发生弹道输运时,当地热力学条件远离平衡态,传统的在热平衡态下定义的温度失去意义[43],而如何在这样的非平衡态下定义温度一直存在争议[44].本文采用等效平衡温度的概念[45,46],认为空间内某个离散单元在发射能量时,产生的声子服从玻-爱因斯坦分布,根据发射能量反解出的分布函数中的温度,便是等效平衡温度.使用DOM离散辐射强度后,等效平衡温度的计算式为其中 σphonon 是声子斯特藩-玻尔兹曼参数[47].使用这种等效平衡温度的好处是,在接近扩散导热时,解声子BTE所得结果能够与傅里叶导热定律的结果相吻合[39].
2.2 拓扑优化
拓扑优化本质上属于大规模0-1整数规划问题,区域内的每个点要么是空洞,要么是实体.这样的离散优化问题在数学上是病态的,解的存在性、收敛性以及求解算法的实现都存在着巨大的困难[48].实现拓扑优化的关键是对优化问题进行适当的正则化处理,将设计变量由离散的变为连续的,从而能够直接使用针对连续设计变量优化问题的求解方法.在这方面,SIMP方法[49]是经典而常用的技巧,其核心思想是将材料性质设定为关于相对密度的幂函数,从而将设计变量转变为连续的相对密度.考虑到(2)式作为描述非傅里叶导热的本构方程,仅有声子消光系数 κ 是与材料类型有关的性质[38],于是本文提出对 κ 进行插值的SIMP方法来实现拓扑优化.空间中任意位置的消光系数写作
其中 ρ 是值在[0,1]区间上连续分布的相对密度变量; κs 和 κf 分别是基底材料和填充材料的消光系数,亦即各自声子平均自由程的倒数; p是值大于1的惩罚因子.惩罚因子的作用是在优化过程中对介于0和1之间的中间密度值进行惩罚,使得 ρ 的值逐渐趋向于0和1两个端点值,从而在优化结果中得到清晰的空洞(ρ=0)和实体(ρ=1)分布.关于SIMP方法插值方式的选择,将在下一节结合优化结果做进一步讨论.
即使采用SIMP方法,拓扑优化仍会面临解的不存在性、不收敛性和不唯一性等问题,导致优化结果会出现棋盘格、网格依赖性等问题[50].为了得到可靠的优化结果,需要引入更多的正则化处理以抑制数值不稳定性.本文采用对相对密度变量 ρ 的全局梯度施加显式约束的方法[15],即在目标函数中增加关于 ρ 的全局梯度的权重惩罚项,以约束ρ的变化程度.至此,可写出通过求解EPRT来获得温度分布的二维体点导热拓扑优化的数学模型为：
其中 Tave 是系统平均温度,Tave,ref 是人为选取的参考温度,α 是相对密度变量全局梯度的权重系数,是无量纲的相对密度变量梯度,A=a2 是二维区域的面积,ϕ 是预先给定的高导热填充材料占整个区域的比例.对平均温度和设计变量梯度进行无量纲化的目的是使二者通过线性求和构成实际的目标函数g.现在的模型以降低系统平均温度为优化目标,如果目标是减小系统的最高温度,则可以将优化对象变更为文献[51]中提到的几何平均温度,其他设置不变.需要说明的是,因为本文的关注点是如何在导热优化中考虑声子弹道输运并分析弹道输运对优化结果的影响,所以优化模型中并未考虑不同材料之间的界面热阻,这样可以更直观地反映弹道效应的影响.界面热阻对现有优化结果的影响将在第3节中进行讨论.
求解(5)式的拓扑优化的流程图如图2所示,包括以下步骤：1)初始化,对区域进行离散并任意给
定一组 ρ 的初始值; 2)系统重分析,根据SIMP方法插值得到当前设计点对应的消光系数,再求解EPRT和稳态能量守恒方程获得声子辐射强度,进而计算平均温度和目标函数; 3)灵敏度分析,计算目标函数和约束函数对设计变量的导数; 4)优化求解,根据导数信息,确定满足约束且减小目标函数的设计变量演化方向(可行下降方向),并计算最佳步长,更新当前设计变量值; 5)收敛判断,满足收敛条件时结束迭代,输出结果; 否则重复步骤2)-4).

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。