, 朱煜, 成荣, 李鑫, 王磊杰, 胡楚雄清华大学 机械工程系, 摩擦学国家重点实验室, 精密超精密制造装备及控制北京市重点实验室, 北京 100084
收稿日期：2021-04-28
基金项目：国家自然科学基金（51677104，51922059）；北京市自然科学基金（JQ19010）
作者简介：赵家琦(1996—), 男, 博士研究生
通讯作者：张鸣, 副研究员, E-mail: Zhangm@tsinghua.edu.cn
摘要：为了提高平面电机水冷散热器热流性能，该文构建了散热器热流模型并开展了"结构尺寸-流道拓扑"并行优化设计。首先，构建了散热器3层热流模型，该模型囊括了盖板与流固混合层的流动/传热耦合效应以及散热器厚度方案对热流场的影响。其次，基于3层热流模型，引入孔隙度域来描述流道拓扑，构建了包含流道拓扑变量与层厚变量的连续伴随结构优化模型，并提出了相应的散热器层厚与流道拓扑并行优化方案。最终，数值案例部分评估了3层热流模型准确性及计算效率、"层厚-流道拓扑"并行优化设计结果热流性能和"尺寸-拓扑"并行优化方案。结果表明：相较于全3D热流模型，所构建的3层热流模型仅需少于10%的计算时间获得了近乎一致的温度场；相较于基准设计，所优化的散热器结构新颖且展现出高达30.82%的目标性能改善；相较于离散优化方案，所提出的并行优化方案高效且取得了具有竞争力的设计。
关键词：散热器热流建模多层模型拓扑优化尺寸优化并行优化
Thermofluid modeling for concurrent size-topology optimization of heat sinks for planar motors
ZHAO Jiaqi, ZHANG Ming
, ZHU Yu, CHENG Rong, LI Xin, WANG Leijie, HU ChuxiongBeijing Key Laboratory of Precision/Ultra-precision Manufacturing Equipments and Control, State Key Laboratory of Tribology, Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China
Abstract: The thermal-hydraulics of water-cooled heat sinks for planar motors were analyzed using a thermofluid model of a heat sink to concurrently optimize the heat sink size and topology. The three-layer thermofluid model of the heat sink included the coupled flow and heat transfer effects between the cover plates and the fluid-solid mixing channel including the influence of the flow channel thickness on the thermal-hydraulics. The model used a porosity model to describe the channel topology and a continuous-adjoint structural optimization model for the geometric variables related to the channel topology and thickness in a concurrent size-topology optimization scheme. Various numerical examples show the accuracy and efficiency of the three-layer thermofluid model and the size-topology optimization scheme. The three-layer model uses 10% less calculational time than full 3D models while still accurately predicting the temperature distribution. The optimized heat sinks have unique topologies with up to 30.82% better heat transfer than baseline designs. This concurrent approach is efficient and obtains designs that are competitive with discrete optimization approach results.
Key words: heat sinkthermofluid modelingmulti-layer modeltopology optimizationsize optimizationconcurrent optimization
随着高性能电子器件朝着集成化、小型化方向不断发展，高功率密度器件的散热技术引起了人们极大的研究兴趣[1]。平面电机广泛应用于以光刻机双工件台为代表的集成电路制造领域[2]。平面电机散热性能直接影响超精密运动系统的加速度与运动精度，从而进一步影响芯片产率与套刻精度[3-4]。现有光刻机工件台中利用平面式水冷散热器带走绝大部分的平面电机线圈发热[4]。结构设计方案对散热器的热流性能具有重要的影响[4-5]。特别地，对流散热器流道拓扑与层厚方案直接影响速度场分布，从而影响流体相与固体相共轭传热，最终影响散热器热耗散性能。此外，流道拓扑与层厚方案直接影响散热器流动耗散性能，从而影响散热器工作泵浦功率。因此，本文工作旨在优化设计平面电机散热器结构(流道拓扑与层厚方案)，以获得更好的目标热流性能。
结构拓扑优化方法在给定的目标与约束下优化设计域内各处材料分布，具有极大的设计自由度，并且不依赖于设计者的经验[6]。自2009年热流问题拓扑优化方案[7]被首次提出，已引起了广泛的研究兴趣。尤其，拓扑优化方法在对流散热器优化设计问题中得到了诸多应用与拓展，包括：层流/湍流、稳态/瞬态、强制对流/自然对流、2D/3D等问题[8-9]。
现有的拓扑优化研究中，热流问题建模方法主要有3种：全3D模型、2D模型与简化模型[8]。全3D模型普适性好，准确度较高，但具有计算量大的特点。特别是在热流问题拓扑优化研究中，几十/上百次的迭代过程使得全3D模型的计算成本极其高昂[10-11]。2D模型计算量小，在早期拓扑优化研究中被广泛采用。但是，2D热流模型没有考虑第三维度的流动与传热效应，计算准确度过低。简化模型的思路是依据散热器结构特点提出适当假定来简化流动与传热模型。其中，多层模型的建模思路被广泛用于平面式对流散热器[10-14]。在该建模思路中，3D热流问题被建模为每一层的2D热流问题以及层间的热流耦合问题。相较于2D模型，多层模型获得了更精确的热流场。相较于3D模型，多层模型具有更低的计算成本。
现有平面式对流散热器多层热流模型研究[10-14]中，相关模型均建立在预设层厚方案下，而没有考虑层厚方案对散热器热流模型的影响。进一步地，现有相关研究[10-14]仅涉及散热器流道拓扑优化而未考虑散热器层厚方案优化。因此，本文旨在开展进一步研究，以优化平面式对流散热器流道拓扑与层厚方案。
基于已有多层热流模型研究，优化层厚方案的基本方法是在有限的层厚方案下优化流道拓扑，并选择相对最佳的设计方案。这种基于穷举的离散方案计算量较大。同时，优化结果严重依赖于设计者预设的层厚方案[15]。相较于离散方案，并行优化方法具有更高的优化效率，能够取得高质量的解[15-16]。然而，现有多层热流模型未考虑层厚设计因素，并行优化方案难以应用于层厚—流道拓扑集成优化设计。
本文构建了包含层厚方案的平面式对流散热器多层热流模型，并介绍了散热器厚度与流道拓扑并行优化方案。首先，构建了平面式散热器的流动/传热3层模型。3层传热模型考虑了盖板传热效应与层厚设计因素。在此基础上，构建了以流道拓扑与散热器层厚为变量、考虑散热性能与流动耗散的多目标优化模型。基于灵敏度信息，构造了流道拓扑变量与层厚变量的并行优化实施方案。数值案例部分展示了所构建的3层热流模型的高效性与准确度，探究了不同目标权重对优化的散热器拓扑的影响。相较于基准设计，所优化的结构取得了高达30.82%的目标性能改善。同时，与离散优化方法相比，并行优化方案高效地获得了具有竞争力的结构拓扑方案。
1 平面式对流散热器热流建模1.1 概述本文所优化的平面式对流散热器如图 1所示。散热器由高导热材料制成。热源(线圈)位于底板下方，冷却介质(水)通过共轭传热效应(对流与传导)增强热源热量耗散。散热器包括3部分：底板、盖板与流固混合区域。散热器层厚分配方案与流固混合区域拓扑方案都对其热流性能有显著的影响。
 |
| 图 1 平面式对流散热器及热源示意图 |
| 图选项 |
散热器入口为流量入口条件，出口压强被定义为0 Pa。所有的流固界面被定义为无滑移边界条件。入口水温为Tin，除热源与入口外，散热器其余外表面被定义为热绝缘。
本文采用多层模型思路建立平面式散热器的流动与传热模型。如图 2所示，底板、流固耦合区域与盖板分别被简化为层1、层2与层3。层1—3均为相应区域的Z向中面。散热器流动问题被降维至层2的流动问题。散热器传热问题被降维至每一层的传热问题以及层间耦合传热问题。
 |
| 图 2 三层热流模型示意图 |
| 图选项 |
1.2 流动建模本文所考虑的流动问题为稳态下黏性不可压缩层流。流体属性被假定为与温度无关。在层2中，流体流动的连续性方程与运动方程分别为
| $\rho_{\mathrm{f}}(\nabla \cdot \boldsymbol{u})=0,$ | (1) |
| $(\boldsymbol{u} \cdot \nabla) \boldsymbol{u}=-\frac{1}{\rho_{\mathrm{f}}} \nabla \boldsymbol{p}+\mu_{\mathrm{f}} \nabla^{2} \boldsymbol{u}-\frac{12 \mu_{\mathrm{f}} \boldsymbol{u}}{d_{2}^{2}}.$ | (2) |
如节1.1所述，入口流量边界条件的控制方程为
| $m=-\int \rho_{\mathrm{f}}\left(\boldsymbol{u} \cdot \boldsymbol{n}_{\mathrm{in}}\right) d_{2} \mathrm{~d} \varGamma_{\mathrm{in}}.$ | (3) |
无滑移边界的流速被定义为0，其控制方程为
| $\boldsymbol{u}=0.$ | (4) |
层1：层2液相下方区域
| $\nabla \cdot\left(d_{1} k_{\mathrm{s}} \nabla \boldsymbol{T}_{1}\right)+Q_{\mathrm{s}}+h_{12}^{\mathrm{f}}\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{1}\right)=0,$ | (5.1) |
| $\nabla \cdot\left(d_{1} k_{\mathrm{s}} \nabla \boldsymbol{T}_{1}\right)+Q_{\mathrm{s}}+h_{12}^{\mathrm{s}}\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{1}\right)=0 ;$ | (5.2) |
| $\begin{gathered}d_{2} \rho_{\mathrm{f}} C_{\mathrm{f}} \boldsymbol{u} \cdot \nabla \boldsymbol{T}_{2}-\nabla \cdot\left(d_{2} k_{\mathrm{f}} \nabla \boldsymbol{T}_{2}\right)= \\h_{12}^{\mathrm{f}}\left(\boldsymbol{T}_{1}-\boldsymbol{T}_{2}\right)+h_{23}^{\mathrm{f}}\left(\boldsymbol{T}_{3}-\boldsymbol{T}_{2}\right),\end{gathered}$ | (5.3) |
| $-\nabla \cdot\left(d_{2} k_{\mathrm{s}} \nabla \boldsymbol{T}_{2}\right)=h_{12}^{\mathrm{s}}\left(\boldsymbol{T}_{1}-\boldsymbol{T}_{2}\right)+h_{23}^{\mathrm{s}}\left(\boldsymbol{T}_{3}-\boldsymbol{T}_{2}\right) ;$ | (5.4) |
| $\nabla \cdot\left(d_{3} k_{\mathrm{s}} \nabla \boldsymbol{T}_{3}\right)+h_{23}^{\mathrm{f}}\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{3}\right)=0,$ | (5.5) |
| $\nabla \cdot\left(d_{3} k_{\mathrm{s}} \nabla \boldsymbol{T}_{3}\right)+h_{23}^{\mathrm{s}}\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{3}\right)=0.$ | (5.6) |
 |
| 图 3 层间传热系数示意图(层2黑色：固相，白色：液相) |
| 图选项 |
相较于2D传热模型，层1—3的传热控制方程中引入了额外的热通量项hjk(Tj- Tk)，k代表层号。该热通量项用于表征散热器厚度方向上的热传递。确定层间传热系数hjk对传热模型准确度有显著影响。文[11]与文[13]提供了详细的层间传热系数推导过程。本文的推导以该工作为参考。
对于层2的流体区域，h12f反映了基板材料与液相的对流换热能力。h12f被定义为[13]
| $h_{12}^{\mathrm{f}}=\frac{\dot{q}_{1}^{\text {wet }}}{A_{1}^{\text {wet }}\left(T_{1}^{\text {wet }}-T_{\mathrm{fl}}\right)} .$ | (6) |
h12s反映了基板翅片基面与翅片间的传导热阻。从翅片底部到翅片中间高度的传导热阻被用来近似h12s[13]。
| $h_{12}^{\mathrm{s}}=\frac{k_{\mathrm{s}}}{0.5 d_{2}}.$ | (7) |
| $h_{23}^{\mathrm{f}}=\frac{\dot{q}_{3}^{\text {wet }}}{A_{3}^{\text {wet }}\left(T_{3}^{\text {wet }}-T_{\mathrm{fl}}\right)},$ | (8) |
| $h_{12}^{\mathrm{s}}=\frac{k_{\mathrm{s}}}{0.5 d_{2}} .$ | (9) |
需要强调的是，式(6)与式(8)求解的是某一层厚分布方案下的层间换热系数。本文考虑了层厚方案对层间换热的影响。基于James等[17]对两平行板间层流传热效应的分析，层间换热系数与流固耦合层的层厚之间被近似为倒数关系。即，已知层厚方案d2*下依据式(6)与式(8)求解的层间换热系数h12f*与h12s*，在层厚方案d2下有
| $h_{12}^{\mathrm{f}}=h_{12}^{\mathrm{f} *} d_{2}^{*} / d_{2},$ | (10) |
| $h_{12}^{\mathrm{s}}=h_{12}^{\mathrm{s} *} d_{2}^{*} / d_{2} .$ | (11) |
| $T_{\text {in }}=T_{0} \text {. }$ | (12) |
| $\boldsymbol{n}_{\varGamma_{e i}} \cdot \nabla \boldsymbol{T}_{i}=0(i=1,2,3) .$ | (13) |
2 尺寸-拓扑并行优化如图 1所示，流道Z方向尺寸小(毫米级)，散热器XY方向尺寸远大于Z方向尺寸。考虑到加工难易程度与计算成本，本文仅考虑流道的XY方向拓扑优化设计。此外，本文考虑层1—3的层厚分布方案优化设计。
2.1 尺寸-拓扑优化模型基于密度法的基本思想，引入孔隙度域γ来描述内部流道拓扑。每一个网格的γ在0—1之间连续变化。γ=1代表液相，γ=0代表固相。节1中的热流模型用于评估确定的结构设计方案。为了评估流道拓扑由γ描述的散热器的热流性能，本节中提出修正的热流控制方程。
层2的修正流动控制方程为
| $\rho_{\mathrm{f}}(\nabla \cdot \boldsymbol{u})=0,$ | (14) |
| $(\boldsymbol{u} \cdot \nabla) \boldsymbol{u}=-\frac{1}{\rho_{\mathrm{f}}} \nabla \boldsymbol{p}+\mu_{\mathrm{f}} \nabla^{2} \boldsymbol{u}-\frac{12 \mu_{\mathrm{f}} \boldsymbol{u}}{d_{2}^{2}}+\alpha(\gamma) \boldsymbol{u}.$ | (15) |
层1—3的修正传热控制方程为
层1：
| $\nabla \cdot\left(d_{1} k_{\mathrm{s}} \nabla \boldsymbol{T}_{1}\right)+Q_{\mathrm{s}}+h_{12}(\gamma)\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{1}\right)=0.$ | (16.1) |
| $\begin{aligned}&\gamma d_{2} \rho_{\mathrm{f}} C_{\mathrm{f}} \boldsymbol{u} \cdot \nabla \boldsymbol{T}_{2}-\nabla \cdot\left(d_{2} k(\gamma) \nabla \boldsymbol{T}_{2}\right)= \\&\ \ h_{12}(\gamma)\left(\boldsymbol{T}_{1}-\boldsymbol{T}_{2}\right)+h_{23}(\gamma)\left(\boldsymbol{T}_{3}-\boldsymbol{T}_{2}\right).\end{aligned}$ | (16.2) |
| $\nabla \cdot\left(d_{3} k_{\mathrm{s}} \nabla \boldsymbol{T}_{3}\right)+h_{23}(\boldsymbol{\gamma})\left(\boldsymbol{T}_{2}-\boldsymbol{T}_{3}\right)=0 .$ | (16.3) |
| $k(\gamma)= \begin{cases}k_{\mathrm{f}}, & \gamma=1, \\k_{\mathrm{s}}, & \gamma=0 ;\end{cases}$ | (17.1) |
| $h_{i j}(\gamma)= \begin{cases}h_{i j}^{\mathrm{f}}, & \gamma=1, \\h_{i j}^{\mathrm{s}}, & \gamma=0 .\end{cases}$ | (17.2) |
| $\alpha(\gamma)=\alpha_{\max }\left(1-\gamma^{c_{1}}\right),$ | (18.1) |
| $k(\gamma)=\frac{\gamma\left(k_{\mathrm{f}}\left(1+c_{2}\right)-k_{\mathrm{s}}\right)+k_{\mathrm{s}}}{1+c_{2} \gamma},$ | (18.2) |
| $h_{i j}(\gamma)=\frac{\gamma\left(h_{i j}^{\mathrm{f}}\left(1+c_{3}\right)-h_{i j}^{\mathrm{s}}\right)+h_{i j}^{\mathrm{s}}}{1+c_{3} \gamma}.$ | (18.3) |
需要强调的是，节1所介绍的热流边界条件在本文优化模型中保持不变。
2.2 优化问题与实施2.2.1 优化问题在流量入口边界条件下，优化目标旨在散热器以较小的流动耗散实现较大的热量耗散。具体地，为了评估散热器热流性能，热源平均温升Trise与出入口压降Δp被当作目标性能。此外，本文引入了设计域ΩD内液相体积比f的约束。
| $T_{\text {rise }}=\int\left(\boldsymbol{T}-T_{\text {initial }}\right) \mathrm{d} \varOmega_{\mathrm{H}} / \int 1 \mathrm{~d} \varOmega_{\mathrm{H}},$ | (19) |
| $\Delta p=\int \boldsymbol{p} \mathrm{d} \varGamma_{\text {in }} / \int 1 \mathrm{~d} \varGamma_{\text {in }}.$ | (20) |
平面散热器尺寸-拓扑优化列式如下：
| $\begin{aligned}&\ {\rm{find}}\ \gamma, \boldsymbol{d}\\&\min \mathrm{obj}=w_{1} \times T_{\text {rise }} / 1 \mathrm{~K}+w_{2} \times \Delta p / 1 \mathrm{~Pa},\end{aligned}$ | (21) |
| ${\rm{s. t.}}\ f=\int \gamma \mathrm{d} \varOmega_{\mathrm{D}} / \int 1 \mathrm{~d} \varOmega_{\mathrm{D}} \leqslant f_{\max },$ |
| $0 \leqslant \gamma \leqslant 1,$ |
| $d_{i \min } \leqslant d_{i} \leqslant d_{i \max }.$ |
2.2.2 并行优化框架为了避免优化框架中的病态数值问题，采用Helmholtz型偏微分方程[18]过滤孔隙度域
| $-r^{2} \nabla^{2} \gamma_{\mathrm{f}}+\boldsymbol{\gamma}_{\mathrm{f}}=\boldsymbol{\gamma} ,$ | (22.1) |
| $-\boldsymbol{n}_{\varGamma_{d}} \cdot\left(-r^{2} \nabla \boldsymbol{\gamma}_{\mathrm{f}}\right)=0.$ | (22.2) |
为了获得清晰的拓扑布局，采用平滑的Heaviside函数[19]投影过滤后的孔隙度域
| $\boldsymbol{\gamma}_{\mathrm{p}}=\frac{\tanh (a b)+\tanh \left(a\left(\boldsymbol{\gamma}_{\mathrm{f}}-b\right)\right)}{\tanh (a b)+\tanh (a(1-b))}.$ | (23) |
层厚变量与拓扑变量的并行优化框架如图 4所示。孔隙度变量经过滤波与投影过程后，修正热流模型基于投影后的孔隙度域和当前层厚变量求解优化目标与优化约束。当未满足收敛准则时，优化算法基于目标函数与约束函数对设计变量的灵敏度信息迭代设计变量，下一次迭代开始。当收敛准则得到满足时，得到优化设计结果。
 |
| 图 4 层厚-流道拓扑并行优化框架 |
| 图选项 |
2.2.3 优化实施基于灵敏度信息的移动渐近法[19]被用于迭代优化变量。链式规则被用来求解目标函数与约束函数对优化变量的灵敏度。以Trise(目标性能之一)对孔隙度变量的灵敏度为例，其推导过程为
| $\frac{\partial T_{\text {rise }}}{\partial \boldsymbol{\gamma}}=\frac{\partial T_{\text {rise }}}{\partial \boldsymbol{\gamma}_{\mathrm{p}}} \frac{\partial \boldsymbol{\gamma}_{\mathrm{p}}}{\partial \boldsymbol{\gamma}_{\mathrm{f}}} \frac{\partial \boldsymbol{\gamma}_{\mathrm{f}}}{\partial \boldsymbol{\gamma}} \text {. }$ | (24) |
为了避免获得较差的局部解，本文引入了参数延续性策略来缓和优化模型的非线性。整个优化过程被分为3个连续的延续步。当前延续步的优化结果作为下一个延续步的初值，直至延续步3结束。在每一个延续步中，插值函数和投影函数的参数是不同的。每一个延续步的收敛准则被定义为KKT(Karush-Kuhn-Tucker)条件下残差向量的范数不超过0.001或迭代次数超过100。
在数值案例部分，采用有限元软件COMSOL Multiphysics求解热流控制方程。整个优化过程在MATLAB中实现。
3 数值案例3.1 参数设定图 5展示了散热器及热源的XY方向几何尺寸。红色矩形内区域为流道拓扑设计域。散热器总层厚为2 mm。优化流固耦合层厚度d2。盖板与底板的厚度被定义为相等。固/液相材料属性如表 1所示。优化问题参数如表 2所示。整个优化过程被分为3个延续步。参考前人研究[10, 14-15]，本文在各个延续步中采用不同的固相反渗透率αmax、插值函数凸性参数c1、c2、c3，以及陡峭度参数a。各个延续步的参数设定如表 3所示。5个参数数值逐渐增大的设置是为了在优化初期确保更凸的优化问题，而后续延续步中逐渐增加对中间密度的惩罚与固液界面的锐度[10]。上述5个参数的具体数值参照已有工作设定。移动渐近法的参数采用默认参数[20]。孔隙度变量和层厚变量的最大步长分别设定为0.1与0.05。
 |
| 图 5 散热器与热源几何尺寸(单位：mm) |
| 图选项 |
表 1 固相与液相材料属性
| 属性 | 数值 | 属性 | 数值 | |
| kf | 0.61 W/(m·K) | kf | 7 W/(m·K) | |
| ρf | 1 000 kg/m3 | Cf | 4 180 J/(kg·K) | |
| μf | 0.001 Pa·s |
表选项
表 2 优化问题参数
| 属性 | 数值 | 属性 | 数值 | |
| m | 5×10-3 kg/s | T0 | 293.15 K | |
| d2min | 0.5 mm | d2max | 1.5 mm | |
| Qs | 1.2×105 W/m2 | r | 0.5 mm | |
| b | 0.5 | Tinitial | 293.15 K |
表选项
表 3 各延续步参数
| 延续步 | αmax/((Pa·s)·m-2) | c1, c2, c3 | a |
| 1 | 1×106 | 3 | 5 |
| 2 | 3×106 | 5 | 10 |
| 3 | 9×106 | 7 | 15 |
表选项
图 6展示了2种具有典型流道拓扑的散热器结构示意图。图 6a中散热器的3D仿真数值解被用来求解h12f*与h23f*。该散热器热流耦合层层厚为1 mm。流道设计域内圆形扰流柱直径为1.31 mm，6×8均布(X方向间距8 mm，Y方向间距3 mm)。在d2*=1 mm时，求得的层间换热系数h12f*=2 922.7 W/(m2·K) 与h23f*=2 141.2 W/(m2·K)。
 |
| 图 6 两种具有典型流道拓扑的散热器 |
| 图选项 |
3.2 三层模型评估基于节3.1的优化设计问题，本节评估所构建的多层热流模型的准确度与计算效率。一方面，流动模型与传热模型的建模误差会累积在所求解的温度场。另一方面，本文尤其关注热源的温度分布。因此，本节以所求解的热源温度场来评估3层热流模型准确度。
图 7展示了3D模型与3层模型所分析的均布扰流柱散热器的热源温度场。可以看出，2种模型所求得的温度场具有较好的一致性。为了进一步评估2个温度场的相似性，引入已有工作[11]中的相似性指标：
| $g_{\mathrm{T}}=\int \frac{\left|\boldsymbol{T}_{3 \mathrm{~L}}-\boldsymbol{T}_{3 \mathrm{D}}\right|}{\boldsymbol{T}_{3 \mathrm{D}}} \mathrm{d} \varOmega_{\mathrm{H}} / \int 1 \mathrm{~d} \varOmega_{\mathrm{H}} \cdot$ | (25) |
 |
| 图 7 不同热流模型求解的热源温度场 |
| 图选项 |
其中：T3L与T3D分别为3层模型与3D模型求得的热源温度场。gT越小，代表两温度场之间的差异越小。依据式(25)，图 7中两温度场相似性指标为gT=0.61%。
为了评估所构建的3层模型在不同层厚方案下的准确度，采用3层模型与3D模型分析不同厚度方案的均布扰流柱散热器结构。表 4列出了不同厚度方案下2种模型求解的热源温度场的相似性指标。
表 4 均布扰流柱结构热源温度场相似性指标
| d2/mm | 0.5 | 0.7 | 0.9 | 1.1 | 1.3 | 1.5 |
| gT/% | 0.39 | 0.76 | 0.53 | 0.40 | 0.46 | 0.70 |
表选项
进一步地，为了评估所构建的3层模型在不同拓扑方案下的准确度，采用3层模型与3D模型分析图 6b中具有平行直翅片拓扑的散热器结构(翅片宽40 mm，高2.2 mm，间距6.5 mm)。表 5列出了不同层厚方案下2种模型求解的热源温度场相似性指标。
表 5 平行直翅片结构热源温度场相似性指标
| d2/mm | 0.5 | 0.7 | 0.9 | 1.1 | 1.3 | 1.5 |
| gT/% | 0.69 | 0.59 | 0.49 | 0.61 | 0.87 | 1.10 |
表选项
由上述分析可知，所构建的3层模型与3D模型所求解的热源温度场在不同拓扑或层厚方案下均展现较好的一致性。此外，本文所有的计算工作都在同一硬件平台(Intel (R) Xeon (TM) E5-2690 CPU、80.0 G RAM)上进行。相较于3D模型，所构建的3层模型平均节省了90%以上的计算时间。综上所述，所构建的3层模型能够高效、有效地评估所研究的平面式对流散热器的热流性能。
3.3 层厚-拓扑并行优化设计基于所构建的结构优化模型，在4种权重方案下并行地优化散热器层厚方案与内部流道拓扑。设计域液相相对体积比上限设定为fmax=0.8。层厚d2的初始值为1 mm。内部流道设计域的初始投影孔隙度为0.8。表 6列出了4种权重方案下所优化的结构的层厚与目标热流性能。以所优化的孔隙度域中孔隙度为0.5的等值线作为固液边界，构建优化结果的几何模型。图 8展示了所优化的内部流道拓扑。
表 6 不同权重方案下所优化的结构层厚与性能
| 权重因子 | d2/mm | Trise/K | Δp/Pa | |
| w1 | w2 | |||
| 1 | 10 | 1.50 | 105.52 | 8.78 |
| 1 | 1 | 1.12 | 80.20 | 20.23 |
| 10 | 1 | 0.63 | 55.76 | 106.50 |
| 30 | 1 | 0.50 | 48.49 | 217.46 |
表选项
 |
| 图 8 不同权重方案下所优化的孔隙度域 |
| 图选项 |
在入口流量固定的边界条件下，减小层厚d2提高了流场平均速度，增强了液相与固相的对流传热，但也增加了维持流动所需的压降。从表 6可以看出，随着Trise的相对权重逐渐增大，所优化的结构的d2与Trise均逐渐减小，Δp逐渐增大。这与上述分析是一致的。
如图 8所示，所优化的孔隙度域展现了清晰的拓扑，所优化的固体特征具有流线形边界。在Δp的权重相对较大时，所优化的结构拓扑简单，流阻较小(图 8a)。在Trise的权重相对较大时，流道拓扑较复杂。复杂的流道拓扑有利于液相流经热源外围区域以获得更好的散热性能，但同时也增大了结构流阻。
为了评估所优化的流道拓扑的性能，以图 6所示2种经典流道拓扑作为基准。在每个权重方案下，经典流道拓扑散热器的层厚方案取值设定为与该权重方案下所优化的层厚相同。表 7列出了不同权重方案下经典流道拓扑散热器与所优化的散热器的目标性能。可以看出，在每种权重方案下，所优化的散热器均比经典流道拓扑散热器表现出更好的目标性能。特别地，在w1=1，w2=10时，拓扑优化散热器的目标性能相较于均布扰流柱散热器的目标性能改善了30.82%。
表 7 不同权重方案下结构目标性能对比
| 权重因子 | 拓扑优化 | 均布扰流柱 | 平行直翅片 | |
| w1 | w2 | |||
| 1 | 10 | 193.3 | 279.4 | 207.7 |
| 1 | 1 | 100.4 | 116.1 | 104.7 |
| 10 | 1 | 664.1 | 748.8 | 708.3 |
| 30 | 1 | 1 672.1 | 1 901.2 | 1 830.3 |
表选项
3.4 并行方案评估为了评估节2.2中的并行优化方案，以权重方案w1=1, w2=1为例，采用离散方案优化层厚d2与内部流道拓扑。具体地，在不同的层厚方案下优化内部流道拓扑。图 9展示了不同层厚方案下优化得到的内部流道拓扑及其目标性能。
 |
| 图 9 不同层厚方案下所优化的拓扑及其目标性能 |
| 图选项 |
相较于固定层厚下流道拓扑优化模型，“层厚-流道拓扑”并行优化模型仅增加一个尺寸变量。因此，离散手段与并行手段的每次迭代耗时近乎一致。然而，在上述案例中，并行优化案例的迭代次数为46次，5个离散案例总迭代次数为286次。由此可见，并行优化方案较离散优化方案高效。
此外，与5个离散优化结果中优化目标最优值对比，并行优化结果的优化目标降低3.46%。在离散优化方案下，为获得更优的解，需要遍历更多的层厚方案并为此付出更大的计算成本。由此可见，本文所提出的并行优化方案高效地取得了具有竞争力的结果。
4 结论本文研究了平面电机对流散热器的热流多层模型与尺寸-拓扑并行优化设计。首先，考虑到高导热材料制成的盖板的传热效应，构建了包括底板、流固耦合层和盖板的3层热流模型。散热器三维热传导问题被建模为每层的2D传热问题以及层间耦合传热问题。特别地，所构建的3层传热模型考虑了层厚方案的影响。进一步引入孔隙度域表征内部流道拓扑，构建了包含层厚变量与流道拓扑变量的优化模型，介绍了相应的并行优化框架与实施方案。
数值案例部分评估了所构建的3层热流模型的准确度与计算成本，展示了不同权重方案下所优化的结构拓扑与目标性能，并评估了所展示的并行优化方案。数值案例显示：相较于3D模型，所构建的3层模型高效且有效地评估了散热器的热流性能；在不同权重方案下，所优化的层厚变化合理；相较于2种具有典型流道拓扑的基准设计，所优化的散热器拓扑具有更好的目标性能；与离散方案相比，所构建的层厚-拓扑并行优化方案高效地取得了具有竞争力的设计。
本文基于层流传热有限元模型评估了多层热流模型的准确性及结构设计方案的目标性能。后续将开展实验验证工作。
参考文献
| [1] | MOORE A L, SHI L. Emerging challenges and materials for thermal management of electronics[J]. Materials Today, 2014, 17(4): 163-174. DOI:10.1016/j.mattod.2014.04.003 |
| [2] | 袁丁. 多自由度无铁芯永磁式平面电机的建模与设计方法研究[D]. 北京: 清华大学, 2019. YUAN D. Modeling and design method for multi-DOF ironless permanent magnet planar motor[D]. Beijing: Tsinghua University, 2019. (in Chinese) |
| [3] | 曹家勇, 朱煜, 汪劲松, 等. 平面电动机设计、控制与应用技术综述[J]. 电工技术学报, 2005, 20(4): 1-8. CAO J Y, ZHU Y, WANG J S, et al. Survey of the state of the art in planar motor technology[J]. Transactions of China Electrotechnical Society, 2005, 20(4): 1-8. DOI:10.3321/j.issn:1000-6753.2005.04.001 (in Chinese) |
| [4] | ZHANG M, XU Q M, CHENG R, et al. Radiator optimization design for planar motors based on parametric components[J]. Journal of Beijing Institute of Technology, 2020, 29(2): 222-231. |
| [5] | AHMED H E, SALMAN B H, KHERBEET A S, et al. Optimization of thermal design of heat sinks: A review[J]. International Journal of Heat and Mass Transfer, 2018, 118: 129-153. DOI:10.1016/j.ijheatmasstransfer.2017.10.099 |
| [6] | WANG L K, MANGLIK R M, SUNDéN B. Plate heat exchangers: Design, applications and performance[M]. Boston: WIT Press, 2007. |
| [7] | DEDE E M. Multiphysics topology optimization of heat transfer and fluid flow systems[C]//Proceedings of the COMSOL Conference. Boston, USA: COMSOL, 2009. |
| [8] | ALEXANDERSEN J, ANDREASEN C S. A review of topology optimisation for fluid-based problems[J]. Fluids, 2020, 5(1): 29. DOI:10.3390/fluids5010029 |
| [9] | DBOUK T. A review about the engineering design of optimal heat transfer systems using topology optimization[J]. Applied Thermal Engineering, 2017, 112: 841-854. DOI:10.1016/j.applthermaleng.2016.10.134 |
| [10] | HAERTEL J H K, ENGELBRECHT K, LAZAROV B S, et al. Topology optimization of a pseudo 3D thermofluid heat sink model[J]. International Journal of Heat and Mass Transfer, 2018, 121: 1073-1088. DOI:10.1016/j.ijheatmasstransfer.2018.01.078 |
| [11] | ZENG T, WANG H, YANG M Z, et al. Topology optimization of heat sinks for instantaneous chip cooling using a transient pseudo-3D thermofluid model[J]. International Journal of Heat and Mass Transfer, 2020, 154: 119681. DOI:10.1016/j.ijheatmasstransfer.2020.119681 |
| [12] | YAN S N, WANG F W, HONG J, et al. Topology optimization of microchannel heat sinks using a two-layer model[J]. International Journal of Heat and Mass Transfer, 2019, 143: 118462. DOI:10.1016/j.ijheatmasstransfer.2019.118462 |
| [13] | ZENG S, KANARGI B, LEE P S. Experimental and numerical investigation of a mini channel forced air heat sink designed by topology optimization[J]. International Journal of Heat and Mass Transfer, 2018, 121: 663-679. DOI:10.1016/j.ijheatmasstransfer.2018.01.039 |
| [14] | ZHAO J Q, ZHANG M, ZHU Y, et al. Topology optimization of planar cooling channels using a three-layer thermofluid model in fully developed laminar flow problems[J]. Structural and Multidisciplinary Optimization, 2021, 63(6): 2789-2809. DOI:10.1007/s00158-021-02842-1 |
| [15] | ZHAO J Q, ZHANG M, ZHU Y, et al. Concurrent optimization of the internal flow channel, inlets, and outlets in forced convection heat sinks[J]. Structural and Multidisciplinary Optimization, 2021, 63(1): 121-136. DOI:10.1007/s00158-020-02670-9 |
| [16] | ZHAO J Q, ZHANG M, ZHU Y, et al. Concurrent optimization of additive manufacturing fabricated lattice structures for natural frequencies[J]. International Journal of Mechanical Sciences, 2019, 163: 105153. DOI:10.1016/j.ijmecsci.2019.105153 |
| [17] | WELTY J R, WICKS C E, WILSON R E, et al. Fundamentals of momentum, heat, and mass transfer[M]. 5th ed. Hoboken: John Wiley & Sons, 2008. |
| [18] | LAZAROV B S, SIGMUND O. Filters in topology optimization based on Helmholtz-type differential equations[J]. International Journal for Numerical Methods in Engineering, 2011, 86(6): 765-781. DOI:10.1002/nme.3072 |
| [19] | WANG F W, LAZAROV B S, SIGMUND O. On projection methods, convergence and robust formulations in topology optimization[J]. Structural and Multidisciplinary Optimization, 2011, 43(6): 767-784. DOI:10.1007/s00158-010-0602-y |
| [20] | SVANBERG K. The method of moving asymptotes-a new method for structural optimization[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359-373. DOI:10.1002/nme.1620240207 |
