微重力条件下列管相变蓄热器的传热特性研究及拓扑优化

高佳徐　任智彬　赵明

摘要：基于相变蓄热器在航天航空领域的应用，对列管式蓄热器建立物理模型和数学模型，采用数值模拟研究其在不同重力及雷诺数下的传热特性，此外引入基于变密度法的拓扑优化方法对蓄热器进行结构优化，设计出一种不等长拓扑肋片模型，并对比研究各因素對相变材料熔化速度的影响。研究结果表明，相同雷诺数时微重力条件下相变材料的熔化速度相比于有重力条件下的较慢，不同雷诺数时微重力条件下的熔化速度比有重力条件下的熔化速度分别减小了 82.7%，86.1%，90.1%。但在微重力条件下，增大雷诺数已不足以有效提升相变材料的熔化速度。与无肋片管相比，经拓扑优化所得肋片模型的强化换热效果明显，在微重力条件下其相变材料熔化时间相比于无肋片管熔化时间缩短了 47.83%，表明所设计的拓扑优化肋片模型可有效地减小微重力对蓄/放热速率的影响。

关键词：微重力；列管式相变蓄热器；数值模拟；拓扑优化

中图分类号：TK 121            文献标志码：A

Heat transfer characteristics and topology optimization of tubularphase change regenerator under microgravity

GAO Jiaxu， REN Zhibin， ZHAO Ming

（School of Energy and Power Engineering， University of Shanghai for Science and Technology， Shanghai 200093， China）

Abstract： Based on the application of phase change thermal storage in the field of aerospace， physicaland mathematic model of the tubular regenerator was built. Numerical simulation was adopted to studythe heat transfer characteristics under different gravity and Reynolds number. A topology optimizationmethod based on variable density method was introduced to optimize the structure of a regenerator. Along range topology fin model was designed， and the influence of various factors on the melting rate ofphase change material （PCM） was compared. The results show that the melting rate of PCM under thesame Reynolds number is slower than that under the gravity condition， and the melting rate under thedifferent  Reynolds  number  is  82.7%，  86.1%  and  90.1%  lower  than  that  under  the  gravity  condition，respectively.  However，  increasing  Reynolds  number  is  not  enough  to  effectively  improve  the  meltingrate  of  PCM  under  microgravity  conditions.  Compared  with  the  unfinned  tube，  the  heat  transferenhancement effect of the fin model obtained by topology optimization is obvious， and the melting timeof the phase change material is reduced by 47.83% compared with the melting time of the unfinnedtubeunder  microgravity  condition，  indicating  that  the  topologically  optimized  fin  model  can  effectivelyreduce the influence of microgravity on the heat storage and release rate.

Keywords： microgravity; tubular phase change regenerator; numerical simulation; topologicaloptimization

航天器在太空轨道运行时，因其外界环境热流变化很大，其内部电子元件设备产生的大量热负荷需要及时排出，以维持设备的稳定运行[1]。若直接利用制冷装置进行废热排放，则势必会使航天器的运行能耗增加，而利用相变蓄热装置既可有效为电子元器件散热，又可減小换热器体积。目前相变热控技术在航空航天领域的应用已有不少成功的案例，如嫦娥一号卫星采用相变材料热管技术，改善了由于月食引起的热环境不稳定问题，极大地降低了温度对光学电子设备的影响[2]。鉴于在微重力条件下，蓄热器中 PCM （相变材料）的熔化过程无法产生自然对流效应，换热流体与相变材料间的整体换热系数下降，从而导致相变蓄热器应用于航天航空领域时效率较低。因此，如何提升微重力条件下蓄热器的蓄热效率成为当前学术界与工程界的重点研究方向。

对于常重力条件下，通常采用的是高热导率的 PCM 和加装各种优化结构的肋片以改善蓄热器结构，达到提高相变蓄热器的整体蓄/放热速率的目的。航天器的结构尺寸受制于运载火箭发射成本和有效载荷运输能力，当蓄热器应用于现代航天器时，其设计尺寸大小和复杂性程度的控制尤为重要[3]。在微重力条件下，添加肋片仍是精简蓄热器整体结构并提升相变蓄热器蓄/放热速率的重要手段。其蓄热器不同的肋片结构应用于微重力条件下的传热特性已有不少学者进行了研究。如甘玉娟等[4]对内嵌矩形和三角形翅片的正十八烷蓄热腔体进行数值模拟，研究了翅片和重力对相变材料熔化速度及传热的影响，结果表明，添加翅片可有效地抵消微重力对相变换热器蓄热效率的影响。麻才新等[5]对内部填充石蜡添加肋片的板式换热器进行数值模拟，研究了不同重力条件下的温度分布和熔化特性，结果表明，微重力增加了相变材料的传热阻力。

拓扑优化是另一种传热优化的重要手段。连续体结构拓扑优化的概念最初应用于解决结构力学问题[6-8]，基本理念是通过在设计域内任意地改变拓扑结构和拓扑联结来实现结构优化。近二十年来，随着其应用领域的扩大，逐渐延伸至传热学领域。近年来，已有不少学者利用拓扑优化方法为相变蓄热器的结构进行优化研究[9-12]。Tian 等[13]对不同相变材料的套管式换热器进行数值模拟，研究了利用不同目标函数和传热机理的仿生拓扑优化来提高蓄热器蓄/放热性能的方法，得出了一种与叶片结构相似的仿生肋片，结果表明，优化后的翅片结构可使熔化时间和凝固时间分别缩短93%和80%。Chen 等[14]以 Ca（OH）2脱水过程为基础对化学蓄热装置中高导电翅片进行拓扑优化，研究了多管热化学蓄热装置中传热不良的问题，拓扑优化后的翅片在水平和垂直方向对称，并且单管的翅片向不同方向延伸，结果表明，拓扑优化的几何结构具有优越的传热性能，与通常使用的纵向设计相比，蓄热速率提高了43%。由于列管式相变蓄热器具有结构紧凑和换热面积大等特点，为拓展其在航空航天微重力领域下的应用，本文采用正十八烷作为相变材料，对列管式相变蓄热器进行数值模拟，研究分析了重力条件及雷诺数对蓄热器蓄热效率的综合影响，并引入拓扑优化方法对蓄热器进行结构优化，经几何重构以筛选出效果最佳的肋片结构。

1物理模型及数学模型

1.1物理模型

列管蓄热单元物理模型如图1所示，模型管长1000 mm，蓄热腔体中列管以正三角形排列，热流体从下方流入，流经装有 PCM 的管束并与 PCM 进行热交换后，从上方流出。各单管为内径40 mm、壁厚1 mm 的圆管，一般认为管间距和管径比值在1.25以上为宜[15]，管间距取60 mm，两者比值为1.5。由于列管排列具有对称性，因此，可将模型简化为二维模型。简化后的二维模型尺寸为60 mm×800 mm ，为确保充分发展的入口速度分布及避免出口回流，拓展进口长度为60 mm ，出口长度为178 mm。相变材料采用正十八烷，其物性参数如表1所示。

1.2数学模型

蓄热器的蓄/放热过程主要是管外流体与管内 PCM 之间的耦合传热，从数值计算角度来说，即为固液相变?导热?强制对流的耦合换热问题。计算区域可分为2部分，一部分为管外强制对流换热流体区，另一部分为管内固液相变区。蓄热壁面采用耦合边界条件。

1.2.1管外换热流体区

流体区域采用湍流模型中的雷诺时均模拟方法，即在时间域上对流场物理量进行雷诺平均化处理，然后求解得到时均化控制方程[16]。本文采用Realizablek-ε湍流模型处理湍流区域。

式中：uf为水流在x方向的速度分量；vf为水流在 y方向的速度分量； P 为压力；ρf为水流密度；cf为水流比热容；uf为水流动力黏度； g为重力加速度；kf为水流导热系数； t为时间；Tf为水温。

1.2.2相变材料区

考虑有重力时的自然对流效应，在 PCM 计算区域同时建立二维层流连续性方程、动量方程和能量方程，相变区域采用焓?多孔模型进行处理。焓法体现在对能量方程的处理上，即直接采用焓作为待求变量。多孔本质上是处理动量方程中的固相?液相?固液界面的统一性问题。液相分数被引入到动量方程的源项中，这样就可将相变界面模糊的“两相区”视为多孔介质区域，并设孔隙率与液相分数相等。当固相孔隙率为零且速度为零时，源项就为无穷大，相当于将固体看成流阻很大的液体。从数学上讲，就是通过源项的作用，将固相区本不应存在的流动速度衰减到最小；对于液相则孔隙率为1，源项为0；对于糊状区，孔隙率则在0～1之间变化，这样由于糊状区凝固时孔隙率会不断地减少，从而引起动量的下降。通过上述源项的作用，就能够将用于液相区的动量方程应用于固相及模糊两相区，从而使方程得到统一，方便数值计算。

式中： u为 PCM 在x方向的速度分量； v为 PCM 在 y方向的速度分量； T 为相变材料的温度； P为压力；ρ为 PCM 密度； c为 PCM 比热容；?为 PCM 动力黏度； k为 PCM 导热系数； t为时间； T0为Boussinesq假设的工作温度；β为 PCM膨胀系数； H为 PCM 的焓。

1.2.4边界条件

传热流体侧：采用给定温度和流速的第一类边界条件[19]，入口温度 Tin=353.15 K ，流体区域初始温度为298.15 K ； PCM 管道壁面采用1 mm 厚的铜材料。流体与壁面交界处采用 coupled 耦合边界条件。

相变材料侧：相变材料区域初始温度为298.15 K ，以保证热流体进入之前相变材料区域为固体。相变材料和肋片交界处采用 coupled 耦合边界条件。

2数值计算方法及模型验证

2.1数值方法

使用 ANSYS Fluent 18.2软件进行流动与传热的数值计算，相变区引入基于焓法的熔化/凝固模型处理相变过程，不同状态下的比热容和导热系数的变化使用函数导入模块（ UDF ）实现。流体区计算采用Realizablek一ε湍流模型。入口边界条件为速度入口，速度范围为0.5～5.0m/s，出口为压力出口，采用 SIMPLE 算法处理压力和速度场的耦合[20]，选择 PRESTO 方法处理压力修正方程， QUICK 算法处理动量方程，二阶迎风算法处理能量方程，上下壁面采用对称边界条件。松弛因子设置如下：压力松弛因子0.5，能量和动量松弛因子分别为0.7和0.9。流动和传热的计算结果输出到 COMSOL 软件进行拓扑优化计算，该拓扑优化模型采用变密度法进行相变材料插值，并利用耦合移动渐进线 MMA 方法进行不同优化目标的稳态求解。

2.2网格验证

为了验证网格数对本文数值模拟的无关性影响，采用网格数为2.5×105，3.4×105，5.2×105和6.4×105的模型进行对比验证，图2（ a ）为网格划分局部示意图。为实现管壁面附近流体流动与传热的准确计算，反映流道内流体真实流动与传热状态，本文在管壁面附近设置了更为细密的边界层网格。图2（b）为不同网格数下温度和液相分数随时间的变化关系。

由图2（ b）可知，在4种不同网格数下的液相分数及 PCM 平均温度变化较小，其中，液相分数计算最大误差为2.50%， PCM 平均温度计算最大误差为0.83%，综合考虑计算精度及计算时间，最终选择3.4×105 W 网格数为本文模型计算网格数。

2.3模型验证

本文根据文献[20]中的实验模型建立了数学模型，依舊采用 FLUENT 18.2软件进行模拟，并检测模型中心点 PCM 的温度变化，模拟方法和松弛因子选择同上。将所得数值模拟结果与实验数据进行对比，如图3所示。

由图3可知，数值模拟的温度分布与文献[21]中的实验数据变化趋势基本一致，与实验相比，计算所得最大误差为1.2%，因此，本文利用该数值方法来研究列管蓄热单元的换热特性是可行的。

3数值模拟结果分析

3.1不同重力条件对蓄放热特性的影响

为了研究重力及雷诺数对列管蓄热单元传热特性的影响，分别在有重力及微重力条件下，Re=27366，136834，273669，对应的入口速度分别为0.5，2.5，5 m/s，研究了速度和重力条件对液相分数及 PCM 平均温度的影响，结果如图4和图5所示。

图4的模拟结果表明，随着雷诺数的增加， PCM 的整体熔化时间减少。根据 PCM 平均温度的变化可知，在同一时刻 Re越大，其对应的平均温度越高，由此可知 Re 的增大，加强了管束外部的湍流扰动，进而提升了外部对流传热特性，从而整体上表现出PCM 熔化速度提升。对比图4与图5可知，在有重力及微重力条件下，两者液相分数及 PCM 平均温度随时间的变化趋势表现一致，但在微重力条件下， PCM 熔化速度远小于有重力条件下的熔化速度，当 Re=27366，136834，273669时，以微重力熔化时间减重力熔化时间，再除以微重力熔化时间的百分率分别为82.7%，86.1%，90.1%。因此，列管外侧流体的流动状态是影响蓄热器充放热效率的重要因素。为了更加细致地研究重力对蓄热单元蓄热效率的影响，本文截取了模型前端部分的流场分布，如图6所示。

如图6的流体区域流场分布所示，在相同雷诺数即相同流速时，流体流经列管绕流后会产生脱体区，该区域内流速较小，并形成回流旋涡。由图6可知，无重力条件下脱体区域的面积比有重力的更大，并且在有重力条件下列管内部熔化的液态 PCM会产生自然对流，从而强化了传热，最终可使有重力时蓄热单元 PCM 熔化速度较微重力时更高。

图7和图8为有重力及微重力条件下不同时刻液相分布图，在有重力时，随着管束内 PCM 开始熔化，液态 PCM 在温差及重力作用下开始产生自然对流，使管内 PCM 熔化过程呈横向不对称；而微重力时管束内没有形成自然对流，因而 PCM 熔化过程是均匀地由外至里。沿流体流动方向靠近入口端的 PCM 熔化最快，这是因为流体从入口流经管束到出口的过程中不断与 PCM 进行换热，导致流体温度沿流向产生一定的温度梯度，越往后端，其与 PCM之间的传热温差越小，因而表现出前端的熔化速度大于后端。

综上所述，不同重力条件和 Re均对蓄热器的蓄放热特性具有重要影响，但在微重力条件下， Re 的影响变得微弱，单纯增大强制换热的力度不足以提升 PCM 的熔化速度，列管蓄热器的结构必须进行进一步优化处理。现将采用拓扑优化方法对该蓄热器模型进行结构优化，以消除微重力的影响。

3.2蓄热器拓扑优化结果分析

优化过程采用 COMSOL 5.4软件进行优化求解，数值优化模型采用移动渐近线方程（MMA ）[22]，优化容差为10?4，单元伪密度初始值为ρi初始值，为0.5，导热材料最大体积分数?=0.05，为了使导热材料边缘更为清晰，取惩罚因子θ=20。优化计算过程中对相变区域蓄热过程平均热流密度进行了统计，结果如图9所示。

由图9可知，随着蓄热过程的进行，热流密度整体呈下降趋势。 PCM 开始熔化时，固态 PCM与管壁之间传热温差较大，因此，导热作用强烈，热流密度较大。随着熔化的进行，管壁面产生的液膜逐渐增厚，导热热阻增大，热流密度也急剧下降[23]。但当液态区域扩展到1600 s 时，热流密度在80 W/m2左右，并维持了较长的时间，因此，拓扑优化时选择相变区域的热源热流均值为80 W/m2。

在拓扑优化过程中，根据上述微重力条件下 PCM 区域温度均匀分布的特点，将单管外壁设置为恒温壁面进行拓扑，因蓄热过程吸收热量，所以，优化区域的热源热流密度设置为?80 W/m2，优化目标为区域内平均温度最大值。拓扑优化结果如图10所示。

如图10（a）模型温度分布所示，采用 MMA 算法进行计算，迭代618次得到最高平均温度为351.89 K。其中，高温区域为拓扑优化的高导热铜材料，并且越靠近铜的区域，温度越高，越靠近中心点的区域，温度越低，由此拓扑所得单管肋片结构如图10（b）所示。

图10（b）为拓扑后模型中铜材料的分布情况，拓扑优化结果为单管内3条长肋片加多条短肋片，但该模型为拓扑理想模型，应用于实际中将过于复杂。蓄热器正趋向于轻量化，优化拓扑后的模型尤为重要。因此，根据拓扑所得长度比例对其进行了简化处理。简化模型如图11所示。

为了对比分析上述拓扑优化所得模型的换热能力，本文还设计了相同总面积的另一组肋片与拓扑模型进行对比，如图12所示。其中，肋片总面积为60 mm2，占 PCM 区域总面积的4.7%。

针对图12所示的3种肋片结构蓄热单元模型，选择入口速度为2.5m/s 进行数值模拟，在微重力条件下，以无肋片结构为基础，研究不同肋片对 PCM 熔化速度的影响，结果如图13和图14所示。

由图13可知，添加肋片后， PCM 的熔化速度得到大幅提升，肋片1，2，3分别使 PCM 熔化时间缩短了47.83%，43.49%和38.3%，这是因为添加肋片后增加了 PCM 与管束壁面的导热面积，从而使蓄热单元的蓄热效率得以提升。图14为不同肋片模型的液相分布。由图14可知，在微重力条件下，没有肋片的模型， PCM 熔化过程是均匀地由管束壁面向中心熔化；而添加肋片后， PCM 的熔化过程是沿着管束壁面及肋片周围开始熔化。与肋片3相比，肋片2中肋片數量较多，因此，在熔化初期，肋片 2 熔化速度较快，但随着熔化的进行，由于肋片 2 较短，PCM 中间出现难熔区域，使熔化时间增加。肋片 3 在刚开始时熔化较慢，但由于其长度较长，因而后期熔化速度增加。由于肋片 1 是长短不同的肋片，其较长的肋片距离 PCM 中心点较近，相比于均匀布置的肋片 2 模型，管内 PCM 初始熔化的径向深度更大，从而导致中心不易熔化的区域温升较快，PCM 熔化速度得到提升[24]。肋片 2 与肋片 3 相比，肋片数量较多，使 PCM 初期熔化速度得到提升，最终可进一步提升 PCM 整体熔化速度。综上所述，采用拓扑优化的肋片 1 可在微重力条件下较好地提高列管式换热器的蓄热效率，可有效弥补微重力条件下蓄热器蓄放热速率。

4结论

通过对列管相变蓄热器进行数值模拟，利用水作为换热介质、正十八烷作为相变材料，研究了不同雷诺数和不同重力条件下蓄热器的蓄热速率变化，并分析得出了不同影响因子对 PCM 熔化速度产生影响的原因。最后引入拓扑优化方法设计了一种不等长拓扑肋片模型，用以改善微重力条件对列管蓄热器效率的影响。综合分析得出结论：

a. 重力条件对蓄热器的蓄放热特性有重要影响。不同 Re 时，微重力条件下的熔化速度比有重力条件下的熔化速度分别减小了 82.7%，86.1%， 90.1%。

b.在微重力条件下，蓄热器内 PCM 区域自然对流作用减弱，且当流体流经管束后其脱体区域增大，因而使得流体与管内 PCM 间的传热恶化。即使增加流体区域的 Re，管后脱体区域的大小仍不会有明显改善，因此，单纯通过增大 Re 的方法并不能有效地提升 PCM 的熔化速度。

c.利用变密度法对蓄热器进行拓扑优化，经几何重构设计出了一种不等长肋片优化模型。该模型在管内增加肋片，弥补了相变材料低导热率的问题，提高了 PCM 管内蓄热速率。计算结果表明，与无肋片管相比，拓扑肋片模型强化换热效果明显，在微重力条件下，添加拓扑肋片后，管内 PCM 熔化时间相比无肋片管缩短了47.83%，因此，采用拓扑优化方法进行肋片的高效设计可有效地弥补微重力对蓄热速率的影响，大幅度地增加列管蓄热器的蓄热速率。

参考文献：

[1]范含林.航天器发展对热控制技术的需求分析[J].航天器工程， 2010， 19（1）：7–11.

[2]邵兴国， 向艳超， 谭沧海.嫦娥一号卫星热控设计中热管的应用及验证[J].航天器工程， 2008， 17（1）：63–67.

[3]马尚君， 刘更， 吴立言， 等.航天器结构的模块化设计方法综述[J].机械科学与技术， 2011， 30（6）：960–967.

[4]甘玉娟， 彭浩， 王一峰.微重力下内嵌翅片蓄热腔体内相变材料的融化传热特性[J].热能动力工程， 2020， 35（12）：108–114.

[5]麻才新， 盛强， 童铁峰.一种空间相变换热器热设计与仿真分析及其改进[J].空间科学学报 ， 2018， 38（3）：409–417.

[6] BENDS?EMP，SIGMUNDO. Materialinterpolation schemes in topology optimization[J]. Archive of Applied Mechanics， 1999， 69（9）：635–654.

[7] BENDS?EMP，SIGMUNDO. Topologyoptimization theory，methods，andapplications[M]. Berlin： Springer，2003.

[8] NEVES M M， RODRIGUES H， GUEDES J M. Optimal designofperiodiclinearelasticmicrostructures[J]. Computers & Structures， 2000， 76（1/3）：421–429.

[9] PIZZOLATO A， SHARMA A， MAUTE K， et al. Topology optimization for heat transfer enhancement in Latent Heat Thermal Energy Storage[J]. International Journal of Heat and Mass Transfer， 2017， 113：875–888.

[10] PIZZOLATO A， SHARMA A， MAUTE K， et al. Multi- scaletopologyoptimizationof multi-materialstructureswith controllable geometric complexity-applications to heat transferproblems[J]. ComputerMethodsinApplied Mechanics and Engineering， 2019， 357：112552.

[11] LUNDGAARD C， ENGELBRECHT K， SIGMUND O. A density-basedtopologyoptimizationmethodologyfor thermalenergystoragesystems[J]. Structuraland Multidisciplinary Optimization， 2019， 60（6）：2189–2204.

[12]游吟， 赵耀， 赵长颖， 等.相变储热单元内肋片结构的拓扑优化[J].科学通报， 2019， 64（11）：1191–1199.

[13] TIANY，LIUXL，XUQ，etal. Bionictopology optimization of finsfor rapid latent heat thermal energy storage[J]. AppliedThermalEngineering， 2021， 194：117104.

[14] CHEN J T， XIA B Q， ZHAO C Y. Topology optimization forheattransferenhancementinthermochemicalheat storage[J]. International Journal of Heat and Mass Transfer， 2020， 154：119785.

[15]韩广顺， 王培伦， 金翼， 等.列管式相变蓄热器性能强化的模拟[J].储能科学与技术， 2015， 4（2）：183–188.

[16] LIU J J， XU C， JU X， et al. Numerical investigation on the heat transfer enhancement of a latent heat thermal energy storagesystemwithbundledtubestructures[J]. Applied Thermal Engineering， 2017， 112：820–831.

[17]田扬， 赵明， 胡明禹，等.加肋旋转对相变蓄热器蓄热性能的影响及场协同分析[J].太阳能学报 ， 2021， 42（3）：395–400.

[18] ZHAO M， TIAN Y， HU M Y， et al. Topology optimization offinsforenergystoragetankwithphasechange material[J]. Numerical Heat Transfer， Part A： Applications，2020， 77（3）：284–301.

[19]张童瑶， 曲明璐， 张娆， 等.相变蓄热器多温位释热特性模拟研究[J].上海理工大学学报， 2020， 42（4）：346–353.

[20]寧利中， 张珂， 宁碧波， 等.侧向加热腔体中的多圈型对流斑图[J].应用数学和力学， 2020， 41（3）：250–259.

[21] ZIVKOVIC B， FUJII I. An analysis of isothermal phase changeof phasechange material within rectangularand cylindrical containers[J]. Solar Energy， 2001， 70（1）：51–61.

[22] SVANBERG K. The method of moving asymptotes-a new method for structural optimization[J]. International Journal forNumericalMethodsinEngineering， 1987， 24（2）：359–373.

[23]崔梦冬， 刘道平， 杨亮， 等.二元有机复合相变材料相变过程热特性实验研究[J].能源研究与信息， 2022， 38（1）：53–60.

[24]李迪， 崔国民， 陈家星， 等.可再生能源多能互补的分布式能源系统两级超结构模型[J].能源研究与信息， 2021， 37（3）：125–133.

（编辑：石瑛）