具有散热通道的传热结构的传热性能与散热通道的材料、形态分布和外部热环境密切相关.当材料和外部热环境一定时,传热结构中的散热通道就对整个结构的热量传递起到主要作用,因此,如何在传热结构中合理布置散热通道对提高结构的传热性能、延长其使用寿命至关重要.在发热量比较集中、空间有限的集成电路中,可采用由高导热材料构成的散热通道直接嵌入到元器件内部,将热量导出到外部环境中,再采用其它高效冷却方式进行冷却,可以有效解决空间限制问题和冷却效率问题.传统的传热结构设计一般基于传热学基本理论及实践经验,然而这种方法很难解决复杂传热结构的设计问题.针对这种复杂模型,可使用结构拓扑优化方法,计算出高自由度的散热通道拓扑形态,在此形态的基础上,通过经验改进以及进一步的尺寸优化获得最佳的结构尺寸.
结构拓扑优化以材料分布为优化对象,在设计空间内寻求最佳的分布方案,以得到结构某种性能最优.在拓扑优化方法中,理论相对成熟且有一定工程应用的方法主要有密度法[1]、均匀化方法[2,3]及Level Set方法[4,5]等.近年来这些方法在传热结构拓扑优化设计中已有所应用.如Bendse等[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中布置高导热材料,形成散热通道,将区域内的热量传送到边界.
 | 图 1 二维热传导问题 Fig. 1 Two-dimensional heat conduction problem |
定义 TTP为结构的散热弱度[16],其中,T和P是离散结构中的节点温度矢量和等效节点载荷矢量,等效节点载荷矢量包括生热率载荷和热流边界引起的等效载荷.结构的散热弱度是结构对外部热载荷的响应,由于散热弱度与热势能具有一致性,当散热弱度最小时,在相同的外部热载荷作用下所能传递的热量最大,使结构具有最佳的温度分布[7].最佳的温度分布体现在结构的平均温度梯度小、平均温度低等方面.因此,可以通过使结构的散热弱度最低来确定热传导问题中高导热材料的最优分布.
在基于密度法的结构拓扑优化设计中,SIMP插值模型是广泛采用的模型,如图 2所示(见下页).在该插值模型中引入了一个伪密度变量的概念,每一个离散单元都对应一个伪密度变量.对于传热结构拓扑优化问题,可使伪密度变量关联的材料属性为离散单元的热传导系数.图 2的SIMP插值关系可以表示为
 | 图 2 SIMP插值曲线 Fig. 2 SIMP interpolation curve |
引入p的目的是在优化过程中避免中间密度单元的形成,使优化后的传热结构中尽可能只存在λmax和λmin两种导热系数的材料.当p<1时,对中间密度具有强化作用;当p>1时,对中间密度具有惩罚作用;当p=1时,不对中间密度进行惩罚或者强化.通常p取2~5范围内的值.
为了便于表达数学模型,引入一个高、低导热材料的热传导系数比h=λmax/λmin.式(1)可以表示为
由单元热传导系数与单元热传导矩阵的线性关系可推出
以结构的散热弱度C=TTP为目标函数,在SIMP插值模型下,其拓扑优化模型可表示为
由上述目标函数、体积约束、平衡方程约束、设计变量约束构造拉格朗日函数,当该构造函数满足KKT条件时,目标函数C具有最小值.综合考虑优化设计变量的边界约束以及热传导矩阵的对称性,对于离散单元应满足KKT条件的表达式为
由(6)式可得
将式(8)和式(9)代入式(7),整理得
将式(10)作为优化设计准则,得出基于该准则下的迭代公式.
当
xi+z)时,xw+1i=(Bwi)ηxwi.
当
=max (xmin,xi-m).
当
min (1,xi+m).
式中,η是为了保证数值计算的稳定性而引入的一个阻尼系数;z为密度变化上限值.
基于最优准则法实现的热传导体结构拓扑优化程序流程如图 3所示.首先建立有限元模型,包括几何模型的建立以及材料、热载荷、边界约束的定义.然后初始化设计变量,将材料参数用设计变量来表示.计算优化模型的目标值以及敏度值.计算满足体积约束条件的Lagrangian体积乘子,并通过当前迭代步中的体积乘子和设计变量来更新下一个迭代步中的设计变量.当设计变量满足收敛条件时,优化过程结束.否则,进入下一个循环过程,直到满足收敛条件.对于迭代过程中的收敛性判定,采用相邻迭代的最大密度变化值为判定依据,表达式为
 | 图 3 优化程序流程图 Fig. 3 Flow of topology optimization |
与以目标变化值为收敛判定条件时相比,采用密度变化值判定条件时,所得到的拓扑形态密度分布更为集中,中间密度单元也相对较少,但迭代次数也相应增多. 1.3 数值不稳定处理方法
在拓扑优化过程中,由于场函数不稳定或数值奇异解的存在导致优化后易出现数值不稳定现象,如棋盘格现象、网格依赖性及局部极值等[17].棋盘格和网格依赖这两种现象一般同时出现在优化结果中,而能够有效去除棋盘格的方法通常也能有效克服网格依赖现象.棋盘格的产生与分析单元的选择有关.研究表明[18,19],合理选择高阶单元或采用非协调元,可有效降低或消除棋盘格.如图 4所示,8节点的高阶单元比4节点的低阶单元在每个边中心点处多出了1个中间节点,使节点温度函数的一阶导数在单元节点更加圆滑,可有效降低棋盘格现象.
 | 图 4 8节点单元 Fig. 4 8 node element |
首先,考虑如图 5所示的设计问题.设计空间是0.1 m×0.1 m的正方形.在内部设计域Ωd存在均匀生热率Q1=6×104 W/m3,在下边界存在Dirichlet边界ΓT,边界温度为T0=0 ℃,边界长度L1=0.01 m,低导热系数λ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 |
其次,考虑在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个工况中都明显出现了高导热材料往高热载区域集中的趋势.
 | 图 7 非均匀生热率载荷Dirichlet边界设计模型Fig. 7 Non-uniform heat generation design model satisfying Dirichlet boundary condition |
|
表 1 生热率载荷 Tab. 1 Uniform heat generation rate (W·m-3) |
 | 图 8 相同体积约束下非均匀生热率载荷Dirichlet 边界设计结果 Fig. 8 Topology forms of non-uniform heat generation models under same volume fraction |
图 9(a)的算例模型是四边Dirichlet边界ΓT,并在内部存在均匀的生热率载荷Qd=3×105 W/m3.由图 9(b)拓扑形态可知,高导热材料从四边向中心区域扩散,并且具有四边对称分布.图 10(a)的算例模型是四角Dirichlet边界ΓT,并在内部存在均匀的生热率载荷Qd=3×105 W/m3.由图 10(b)拓扑形态可知,高导热材料从四角向中心区域扩散,且在模型中出现了十字型的空白区域.图 11(a)的算例模型是四角Dirichlet边界ΓT,并在中心点存在集中生热率载荷Q1=3×106 W/m3.由图 11(b)拓扑形态可知,高导热材料由四角和中心点形成了一个X型的材料分布.图 12(a)的算例模型是四边Dirichlet边界ΓT,并在中心点存在集中生热率载荷Q1=3×106 W/m3.由图 12(b)拓扑形态可知,高导热材料由四角和中心点形成了一个十字型的材料分布.综上所述,图 9(b)和图 10(b)均明显出现了有分歧状的拓扑结构,而图 11(b)和图 12(b)出现了清晰的 X和十型的拓扑结构,这是由于集中生热时,热量流动不易受到分支结构的引导,所以,拓扑结构中高导热材料都集中在生热点与热量出口的最短路径上.
 | 图 9 四边温度边界-均匀生热 Fig. 9 Four sides temperature boundary with uniform heat generation |
 | 图 10 四角温度边界-均匀生热 Fig. 10 Four corners temperature boundary with uniform heat generation |
 | 图 11 四角温度边界-中心生热 Fig. 11 Four corners temperature boundary with centre point heat generation |
 | 图 12 四边温度边界-中心生热 Fig. 12 Four sides temperature boundary with centre point heat generation |
对于薄壁管状体,由于管壁相对于管径而言小很多,所以,此类问题可以采用壳单元进行计算.图 13是薄壁管状体算例模型,在管的表面有尺寸0.04 m×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 |
 | 图 14 薄壁管状体结构拓扑结果(不同视角) Fig. 14 Different views of the topology form |
热传导问题中存在Dirichlet边界和Neumann边界的混合边界的算例模型如图 15所示(见下页).设计空间规格为0.1 m×0.1 m的正方形,左下角有宽度为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.01 m增加到0.10 m时,高导热材料越来越多地集中在热量的入口边界上,如图 16(c)所示. 由此可知,当模型中只有Dirichlet边界和Neumann边界时,高导热材料主要分布在Dirichlet边界和Neumann边界的最短路径上.
 | 图 15 混合边界模型 Fig. 15 Model satisfying both Dirichlet and neumann boundary conditions |
|
|
表 2 工况参数 Table. 2 Parameters of working condition |
 | 图 16 混合边界下的优化结果 Fig. 16 Results of different length of boundary |
应用密度法设计传热结构中散热通道的分布,建立了在稳态热传导下的拓扑优化数学模型,推出了满足KKT条件的传热结构拓扑优化的最优准则法迭代公式,通过对Dirichlet边界和Neumann边界的不同组合下的数值算列,验证了算法的有效性.在拓扑结果的形态上,Dirichlet边界和非集中生热率载荷时,高导热材料分布形态为类似植物分枝结构的拓扑形态;Dirichlet边界和集中生热率载荷时,高导热材料分布在集中载荷点与Dirichlet边界的最短路径上;Dirichlet边界和Neumann边界的混合边界时,高导热材料分布在两种边界的最短路径上,且集中在Neumann边界附近.
首页-富联娱乐-富联中国加盟站
