基于双曲正弦函数的拓扑优化在热传导结构设计中的应用

2020-06-27 04:35薛红军陶才勇党思娜

航空工程进展 2020年3期

薛红军，陶才勇，党思娜

(西北工业大学航空学院，西安 710072)

0 引言

飞机电子设备不断增多，其散热结构是航空工程结构设计过程中必须考虑的重要问题之一[1]，飞机中的传热结构设计既要满足一定的强度和刚度要求以达到特定的功能，又要满足散热要求，目的是保证关键设备不会因为温度过高而出现故障。通过热力学计算设计的传统散热装置，已无法满足对热防护结构具有高性能要求的实际工程需求[2]，因此，寻求一种新的散热装置设计方法成为解决热传导问题的关键。

针对热传导问题，H.Rodrigues等[3-4]和C.Jog[5]将拓扑优化应用在热弹性结构中，初次实现了拓扑优化在热传导中的应用；A.Bejan[6]在热传导的背景下介绍了热传导通道构造理论，建立了所谓的“体-点(VP)问题”或“通道问题”，并讨论了在发热域(例如CPU)中布置固定数量的材料的需要；Li Q等[7]明确地将ESO用于热传导问题。

现阶段应用于热传导结构的拓扑优化方法主要有：渐进结构优化方法(ESO)、双向渐进结构优化方法(BESO)[8]、变密度方法[9]、水平集方法[10]、拓扑导数方法[11]。就目前的研究成果而言，应用比较广泛的拓扑优化方法是变密度法。

基于变密度方法的散热结构的拓扑优化设计主要使用：带惩罚的各向同性固体材料模型[12](Solid Isotropic Material with Penalization，简称SIMP)和材料属性的理性近似模型[13](Rational Approximation of Material Properties，简称RAMP)。F.H.Burger等[14]利用SIMP和移动渐近线(MMAs)的方法，探讨了体-点问题的3D解决方案；T.Van Oevelen等[15]采用SIMP模型演示了用双层简化模型表示全三维解的共轭传热问题的解；G.Marck等[16]使用SIMP方法进行了多目标优化；Li Jiachun等[17]使用了材料特性的理性近似模型(RAMP)，进行了热传导结构优化设计的探索。两种模型对于导热结构的拓扑优化设计都能得到比较好的布局，但是对于“体-点”问题[18]，SIMP方法会产生很多树状末梢结构，这样的结果有一定的缺点，如可制造性差，并且结果并不是最靠近最优值的结果[19]。因此有学者通过RAMP方法得到了一定的改进，能过获取相对更加良好的拓扑结果[20]。但是，这两种方法有其局限性，在解决二维、三维均匀热源的散热问题时，SIMP方法不够精确，RAMP方法迭代步长过多。

本文在变密度法的基础上，建立基于sinh函数的插值模型，与SIMP和RAMP模型的计算结果进行对比，并将此模型应用于机载LRM模块导热拓扑优化设计中进行验证。

1 基于sinh函数插值模型的拓扑优化在热传导结构设计中的模型建立

基于变密度法的拓扑优化的关键是选择合适的插值函数和惩罚技术，将问题的物理量表示为连续设计变量的函数。本文采用基于sinh函数的插值模型[19]进行求解计算。

1.1 基于sinh函数的插值模型

基于sinh函数的插值模型，其实质是在插值模型中应用了双曲正弦函数(sinh函数)：

(1)

可以简化为

(2)

sinh函数插值模型采用双曲正弦函数来对中间密度进行惩罚，其函数关系如图1所示。

图1 sinh插值模型

sinh函数插值模型的微分形式为

(3)

当ρ=0时，

当ρ=1时，

1.2 sinh插值模型灵敏度分析

计算导数最好的办法是使用伴随方法，在此也采用伴随方法推导sinh函数插值模型的目标函数的灵敏度。

本文以热柔度为目标函数，基于sinh函数插值模型的目标函数为

(4)

式中：K0为单元导热系数矩阵。

对于上面给出的热柔度目标函数，通过添加零函数，将其重写为

C=QTT-λT(KT-Q)

(5)

式中：λ为拉格朗日乘数，是一个任意的、但是固定的实数向量。

然后将其对设计变量求导得到：

(6)

当λ满足伴随方程

QT-λTK=0

(7)

目标函数灵敏度可写为

(8)

根据伴随方程，可以直接得到λ=T，因此目标函数灵敏度为

(9)

其中，对于传热系数矩阵对设计变量的导数为

(10)

因此，当传热系数矩阵采用sinh函数插值模型对其进行插值，则目标函数的灵敏度为

(11)

式中：K0为具有单元热传导矩阵，当材料的导热系数为k0时，对于二维的四节点矩形单元，单元的尺寸为单位长度，K0可写为

(12)

对于三维八节点正六面体单元，单元的尺寸为单位长度，单元导热系数矩阵K0可写为

(13)

当惩罚因子p=3时，灵敏度随着设计变量的变化曲线如图2所示，可以看出：所有的灵敏度值

都是负值，也就意味着随着导热材料的增加，目标函数在不断地减小，导热性能更好。

图2 灵敏度变化曲线

为了避免棋盘格现象和网格依赖现象，在计算中采用灵敏度过滤技术，将灵敏度改写为

(14)

式中：Hei为单元i的一个集合；γ为一个比较小的常数，为避免分母为0，可取γ=10-3；单元i中心到单元e的距离Δ(e,i)小于最小滤波半径rmin。

1.3 基于sinh函数插值模型的优化准则

本文以结构的散热弱度最小为目标函数，体积分数、结构平衡方程和设计变量约束为约束条件，建立基于sinh函数插值模型的热传导结构的拓扑优化数学模型为

min:

(15)

式中：C为散热弱度；K为结构的总体传热系数矩阵；F为温度载荷向量；T为节点温度向量；Ve为单元体积；V0为结构初始状态下总的体积；f为体积分数，是优化后的结构总体积与初始状态下总体积之比；ximin为设计变量取值的下限，通常取一个较小的数，如0.001；ximax为设计变量取值的上限，取为ximax=1。

利用拉格朗日乘数法将多为高次多约束问题转换为无约束最优化问题，对应于拓扑优化问题的拉格朗日函数为

(16)

式中：λ1为拉格朗日乘子；λ2，λ3，λ4为拉格朗日乘子向量；bi和ci为松弛因子；x为xi组成的列向量。

(17)

(18)

考虑上式等于0的情况，并且代入C=TTK(X)T有：

(19)

利用导热系数矩阵的对称性得出：

(20)

式中：λ2为列矢量，其取值无限制。取λ2=-2T代入式(20)得出：

(21)

进而得出：

(22)

因此可以得到基于优化准则法的迭代公式为

(23)

式中：m为移动极限，通常取一个较小的正数；η为阻尼系数，引入η的目的是为了保证数值计算的稳定性和收敛性。

(24)

(25)

将式(21)代入式(25)得到：

(26)

通过求解上述方程便可以求解出λ1，通常采用二分法进行求解。

设计变量通过优化准则更新，其收敛性的检查方法如式(27)所示：

|max (xk+1-xk)|<ε

(27)

式中:xk+1为更新的设计变量；xk为旧的设计变量；ε为可行性误差，即每个单元更新后的密度值与前一次的密度值之差的绝对值不超过规定的误差值。

2 模型算法实现以及算例分析

2.1 算法实现

本文在MATLAB上进行程序编译，实现模型功能，整体的实现步骤如图3所示。

图3 基于sinh的热传导结构拓扑优化方法算法流程

2.2 算例分析

2.2.1 二维算例

(1) 体-点问题

算例1：整块板均匀加热，左边界中心处给定温度T=0 ℃，其余边界绝热，其结构示意图如图4所示。

图4 导热板均匀加热，左侧中心处冷却模型示意图

高热导率材料的导热系数为k0=1 W/m·℃的假想材料，低热导率材料的导热系数为kmin=0.001 W/m·℃的假想材料，模型离散网格为180×180。优化参数设置为：体积分数0.4；惩罚因子p=3；采用敏度滤波，滤波半径r=2。分别采用SIMP模型，RAMP模型和基于sinh函数模型对其进行优化计算，其优化结果如图5所示，对比结果如表1所示。

(a) SIMP (b) RAMP (c) sinh模型

图5 导热板均匀加热，左侧中心处冷却模型优化结果

Fig.5 Optimization results of heat-conducting plate evenly heating and cooling at the left center

表1 导热板均匀加热，左侧中心处冷却模型优化数据对比

(2) 均匀热源问题其他边界条件

算例2：整块板均匀加热，左边界给定温度T=0 ℃，其余边界绝热，其结构示意图如图6所示。高热导率材料的导热系数为k0=1 W/m·℃的假想材料，低热导率材料的导热系数为kmin=0.001 W/m·℃的假想材料，优化参数设置为：体积分数0.5；惩罚因子p=3；采用敏度滤波，滤波半径r=3。分别采用SIMP模型，RAMP模型和基于sinh函数模型对其进行优化计算，其优化结果如图7所示，对比结果如表2所示。

图6 导热板均匀加热，左侧边界冷却模型示意图

(a) SIMP (b) RAMP (c) sinh模型

图7 导热板均匀加热，左侧边界冷却模型优化结果

Fig.7 Optimization results of heat-conducting plate evenly heating, and the left edge cooling

表2 导热板均匀加热，左侧边界冷却模型优化数据对比

(3) 集中热源问题优化分析

算例3：平板中心加热，四个角点给定温度T=0 ℃，其结构示意图如图8所示,其优化参数与算例2一致。分别采用SIMP模型，RAMP模型和基于sinh函数模型进行优化计算，优化结果如图9所示，对比结果如表3所示。

图8 导热板的中心加热，四个顶点冷却模型示意图

(a) SIMP (b) RAMP (c) sinh模型

图9 导热板的中心加热，四个顶点冷却模型优化结果

Fig.9 Optimization Results of heat-conducting plate center heated and four apexes cooled

表3 导热板的中心加热，四个顶点冷却模型优化结果数据对比

算例4：板上有5个集中热源，四周边界给定温度T=0 ℃，其结构示意图10如图所示,其优化参数与算例2一致。高热导率材料的导热系数为k0=1 W/m·℃的假想材料，低热导率材料的导热系数为kmin=0.001 W/m·℃的假想材料。优化参数设置为：体积分数0.5；惩罚因子p=3；采用敏度滤波，滤波半径r=2.5。分别采用SIMP模型，RAMP模型和基于sinh函数模型进行优化计算，优化结果如图11所示，对比结果如表4所示。

图10 导热板上5个集中热源加热，四周边界散热模型示意图

(a) SIMP (b) RAMP (c) sinh模型

图11 导热板上5个集中热源加热，四周边界散热模型优化结果

Fig.11 Optimization results of the heating of 5 centralized heat sources on the thermal plate and the heat dissipation at the surrounding boundary

表4 导热板上5个集中热源加热，四周边界散热模型优化结果数据对比

2.2.2 三维算例

(1) 均匀热源问题

算例5：设计域下面均匀加热，上平面四个顶点温度固定T=0 ℃，模型示意图如图12所示。高热导率材料的导热系数为k0=1 W/m·℃的假想材料，低热导率材料的导热系数为kmin=0.001 W/m·℃的假想材料。优化参数设置为：体积分数0.3；惩罚因子p=3；采用密度滤波，滤波半径r=1.4。分别采用SIMP模型，RAMP模型和基于sinh函数模型进行优化计算，优化结果如图13所示，对比结果如表5所示。

图12 下表面均匀加热，上表面四个顶点冷却模型示意图

(a) SIMP (b) RAMP (c) sinh模型

图13 下表面均匀加热，上表面四个顶点冷却模型优化结果

Fig.13 Optimization results of uniform heating of the lower surface and cooling of the four vertices of the upper surface

表5 下表面均匀加热，上表面四个顶点冷却模型优化结果数据对比

(2) 集中热源问题

算例6:正方体中心有集中热源，八个顶点温度固定为T=0 ℃，模型示意图如图 14所示。高热导率材料的导热系数为k0=1 W/m·℃的假想材料，低热导率材料的导热系数为kmin=0.001 W/m·℃的假想材料。优化参数设置为：体积分数0.3；惩罚因子p=3；采用密度滤波，滤波半径r=1.4。分别采用SIMP模型，RAMP模型和基于sinh函数模型进行优化计算，优化结果如图15所示，对比结果如表6所示。

图14 正方体中心加热，八个顶点冷却模型示意图

(a) SIMP (b) RAMP (c) sinh模型

图15 正方体中心加热，八个顶点冷却模型优化结果

Fig.15 Optimization results of cube center heating and eight vertices cooling

表6 正方体中心加热，八个顶点冷却模型优化结果数据对比

通过算例将SIMP模型和RAMP模型与基于sinh函数的插值模型进行对比，得出基于sinh函数的插值模型较SIMP插值模型精确，较RAMP插值模型的迭代步长少，计算效率高。

3 机载LRM模块导热拓扑优化设计

机载LRM(现场可更换模块)模块已成为飞机电子系统设备的硬件标准,其优点是标准程度高、模块功能化强、维修性好、成本低。由于机载LRM模块是由超大规模集成电路构成的，包含大量的VLSI和VHSIC芯片，而这些设备在工作时会产生大量热量，使得LRM模块温度急剧上升，因此，机载LRM模块轻量化的冷却成为设计的首要问题。

由于LRM模块空间的限制，传导冷却方法被优先考虑。此方法是通过将热量从LRM模块芯片传至导热体，再经由导热体传至机架冷板，最后通过通风冷却或通液冷却等散热形式将热量带走。简化后的LRM模块导热模型如图16所示，导热体材料为铜，导热系数为377 W/m·℃。