一种基于瞬态自然对流模型的潜热储能装置高效三维拓扑优化两阶段方法

《International Journal of Heat and Mass Transfer》：A two-stage method for efficient three-dimensional topology optimization of latent energy storage devices with a transient natural convection model

　　

摘要

潜热储热（LHTES）装置利用相变材料（PCM）已在各种应用中得到了广泛使用。为了增强传热，拓扑优化（TO）已被用于生成高导热材料（HCM）的新型结构。然而，由于计算成本高，三维LHTES优化的拓扑优化仍然有限。本研究利用界面过滤结构优化（IFSO）方法保持尖锐边界的特征，跳过早期优化迭代，这些迭代由于设计中的显著灰色性而对性能改善贡献甚微。提出了一种两阶段方法。在第一阶段中，在一个拓扑优化框架内使用导热模型。利用导热优化后的设计，在第二阶段中，将IFSO方法与自然对流模型耦合以进一步改进设计。由于使用了导热模型，第一阶段的计算成本大幅降低，而清晰的初始设计减少了第二阶段所需的迭代次数。因此，三维几何形状的优化变得可行。两阶段优化应用于三个示例：二维矩形域、三维立方域和三维环形域。结果表明，导热优化后的设计优于传统设计。在所有示例中，与导热优化设计相比，第二阶段实现了高达15.0%的温度改善。在所有示例中，自然对流的影响是一致的，它在上部区域加速熔化，而在下部区域减缓熔化。为解决这种不平衡，减少了受影响区域的HCM结构，从而提高了整体性能。

1. 引言

潜热储热（LHTES）装置已广泛应用于多个工程领域，包括聚光太阳能发电（CSP）电站[1]、工业废热回收[2]、供暖通风和空调（HVAC）系统[3]、电池热管理系统（BTMS）[4]和集成电路冷却[5]。LHTES系统中的相变材料（PCM）吸收外部热能并在熔化过程中以潜热形式储存。与显热储能系统相比，相变材料由于潜热的存在提供了更高的有效热储存容量。尽管相变材料具有高储能容量，但其固有的低导热率限制了能量系统的充放电速率。各种方法已被研究以提高整体导热率，包括纳米颗粒[6]、泡沫结构[7]和高导热材料（HCM）结构[8]。在上述方法中，插入HCM翅片是最广泛研究和实施的方法。

许多新型设计，如叶脉[9]、螺旋[10]和分形结构[11]，已被开发为HCM结构以优化性能。近年来，拓扑优化（TO）[12]越来越被应用于HCM结构的设计。以下列出了关于瞬态PCM熔化问题的拓扑优化代表性研究。Pizzolato等人首次提出了瞬态对流熔化问题的拓扑优化框架[13]。Chen等人利用拓扑优化为球形几何形状生成HCM结构[14]。Zhang等人将拓扑优化后的翅片插入壳管式构型中以增强传热[15]。Wang和Qian运用拓扑优化技术研究了慢速和快速放电条件下的优化结构[16]。Sazvar等人将HCM结构优化与体积最小化相结合以实现紧凑型LHTES装置[17]。由于其巨大的设计空间和适应自由形状及拓扑变化的能力，拓扑优化在优化LHTES系统方面已展现出显著的成功。

包含瞬态自然对流模型的LHTES系统的拓扑优化在二维中已被广泛研究，并取得了相当大的成功[13]，[15]，[18]，[19]，[20]。然而，扩展到三维几何形状受到不可行的计算负担的限制。据作者所知，迄今为止只有六项研究调查了三维LHTES问题的拓扑优化，如表1所示。表1中列出的先前文献按拓扑优化过程中使用的模型进行分类。对当前文献的回顾表明，尚未有三维拓扑优化研究在LHTES问题中采用瞬态、相变和对流的物理模型，这可归因于高昂的计算成本。对流模型由速度、压力和温度的耦合控制方程组成。求解这些方程需要迭代方法，导致与导热模型相比显著更高的计算成本。在瞬态问题中，计算成本进一步增加，因为每个时间步都需要迭代求解。在三维问题中，由于空间域中的自由度（DOF）急剧增加，计算成本变得令人望而却步。因此，为将计算成本限制在实际水平内，它们要么将传热模型简化为仅导热[21]，[22]，[23]，要么忽略时间依赖性[18]，[19]，[22]，[24]。然而，省略任何一项都会降低预测精度。如[19]中所报告并在本研究中进一步考察，自然对流在上下区域之间的非均匀熔化速率中起着重要作用。在优化过程中忽略这种效果会导致效率较低的设计。如[16]中所报告，放电时长的变化显著影响优化结构，突显了在稳态模型中省略瞬态效果可能导致重要信息的丢失。此外，在某些特定场景中，例如对绝热LHTES装置施加恒定热通量时，从基本原理上来说不可能达到稳态。因此，考虑到相变的瞬态对流模型对于三维LHTES拓扑优化是必不可少的。

为了有效地对具有全物理模型（瞬态、相变和对流）的三维LHTES问题执行拓扑优化，本研究中提出了一种两阶段优化方法以减少计算时间。在标准的基于密度的拓扑优化过程中，前100-200次迭代生成的设计通常表现出显著的灰色性，因为尖锐度因子在整个优化过程中逐渐增加。然而，灰色设计的速度和温度预测缺乏物理意义，从非物理场计算出的灵敏度可能将优化引向较差的局部最小值。这表明瞬态对流拓扑优化中的早期迭代计算成本高昂，但由于高灰色性对性能改善贡献甚微。如果优化从清晰的初始设计开始，可以避免前100-200次以高灰色性为特征的迭代，从而减少不必要的计算成本。为实现这一目标，界面过滤结构优化（IFSO）方法[25]，一种基于密度的拓扑优化方法的变体，在全物理优化中被采用，因为它在整个优化过程中保持尖锐的设计边界，因此消除了初始高灰色性迭代。

表1. 与先前三维LHTES问题的拓扑优化研究的比较。

文献 | 时间 | 相变 | 模型
Pizzolato等人 2017 [21] | 瞬态 | 凝固 | 导热
Iradukunda等人 2020 [24] | 稳态 | 单相 | 对流
Lum等人 2024 [22] | 稳态 | 单相 | 导热
Guibert等人 2024 [23] | 瞬态 | 熔化 | 导热
Zhang等人 2025 [18] | 稳态 | 熔化 | 对流
Yin等人 2026 [19] | 稳态 | 熔化 | 对流
本文工作 | 瞬态 | 熔化 | 对流

由于不依赖灰色性，IFSO方法需要事先确定一个合适的初始设计。为了获得IFSO过程的初始设计，在本研究中，优化分两个阶段进行：第一阶段采用低保真度模型生成合理的初始设计，然后用作后续全物理模型IFSO阶段的起点。更具体地说，优化流程的流程图如图1(a)所示。两阶段优化框架将整体过程分为导热阶段（第一阶段）和对流阶段（第二阶段）。在导热阶段中，在灰色初始设计上进行基于导热模型的拓扑优化。第一阶段优化后的设计，具有尖锐边界，作为第二阶段的初始设计。将通过数值示例证明，导热优化后的设计优于传统设计，意味着它是一个合适的初始设计。在第二阶段中，采用包含自然对流的全模型在IFSO方法中进一步提升从导热优化设计获得的性能。如图1(b)所示，IFSO方法使用具有非均匀半径的Helmholtz型偏微分方程（PDE）过滤器来诱导边界运动。不同的半径模式（图1(b2)中的1-2和3-4）诱导不同的运动行为（1-2：翅片侵蚀；3-4：翅片膨胀，如图1(b4)所示）。在IFSO方法中，PDE过滤器半径由优化器确定并更新以最小化目标值。设计边界在整个优化过程中保持为尖锐的，如图1(b3)中的示例所示。

研究三个数值示例以说明两阶段方法，并从优化设计中提取自然对流效应的物理见解。两阶段方法首先应用于二维矩形几何形状，其中第一阶段导热优化设计的性能比传统设计高出2.4%。在第二阶段中，对流优化进一步将IFSO对流设计的热性能提高了高达15.0%。此外，还研究了两个三维几何形状：立方域和环形域。对流优化阶段实现了与导热优化设计相比平均基底温度降低高达13.3%。在所有示例中，第二阶段观察到一致的形状演化机制：自然对流诱导浮力驱动流动，将热量从底部带至顶部，增强上部区域的熔化速率，同时降低下部区域的熔化速率。为解决熔化速率的不平衡，减少或去除受影响区域的HCM结构，提高整体性能。该机制为自然对流条件下的LHTES装置提供了通用的设计指导。

下载：下载高分辨率图像（313KB）下载：下载全尺寸图像
图1. 内嵌IFSO的两阶段优化策略介绍。

本文其余部分组织如下。第2节介绍了四个不同保真度水平的物理模型。第3节详细介绍了提出的两阶段方法。第4节展示了本研究中考虑的三个数值示例，包括一个二维案例和两个三维案例。最后，第5节以关键发现总结了本文。

2. 物理模型

本节介绍了本研究所研究的潜热储能问题的几何构型、物理模型和边界条件。本研究中使用了四个物理模型，每个模型反映了不同的物理细节和简化程度：完整的自然对流模型（第2.1节）；用于优化的简化自然对流模型（第2.2节）；仅导热模型（第2.3节）；和用于优化的简化仅导热模型（第2.4节）。完整模型与优化中使用的模型之间的比较在2.5节中通过一个数值示例给出。

下载：下载高分辨率图像（242KB）下载：下载全尺寸图像
图2. 优化和CFD模拟的几何构型。

2.1. 完整自然对流模型

用于优化后设计数值验证的完整自然对流模型描述如下。如图2(a)所示，在该模型中，整个域被分为两个子域：PCM域Ω_PCM和HCM域Ω_HCM，由一个明确定义的界面Γ_c分隔。

在HCM域中，能量储存和传输通过导热方程(1)来建模，

(1) ?h_h/?t = ?·(k_h?T_h)

其中h_h是HCM域中的焓，由方程(2)定义

(2) h_h = ∫_T_ref^T ρ_h C_{p,h} dT?。

ρ_h、C_{p,h}和k_h分别是HCM的密度、比热和导热率。

在PCM域内，采用瞬态不可压缩Navier-Stokes方程与瞬态能量方程耦合来预测PCM的流动运动和能量储存。采用Boussinesq近似来简化动量方程中的自然对流项[26]。为考虑固-液相变，采用焓-多孔介质法[27]。在动量方程中，引入Kozeny-Carman关系项来表示PCM的熔化并在固体区域惩罚速度[28]。控制方程列于方程(3)中

(3a) ?·u_p = 0，
(3b) ?u_p/?t + u_p·?u_p = -?p_p + ?·(ν_p u_p) - β_p(T_p - T_ref)g - F(T)u_p，
(3c) ?h_p/?t + u·?h_p = ?·(k_p?T_p)，

其中ν_p、β_p和k_p分别是PCM的运动粘度、热膨胀率和导热率。g是重力加速度，在本研究中等于9.8 m/s2。u_p、p_p和T_p分别是PCM域中的速度、压力和温度。h_p是PCM的焓，由方程(4)定义

(4) h_p = ∫_T_ref^T ρ_p C_{p,p} dT? + ρ_p L_p f(T)，

其中ρ_p、C_{p,p}和L_p分别是PCM的密度、比热和潜热。

区分固相和液相PCM相的液相分数函数??(??)由公式(5)定义：

(5) ??(??) = { 0, ?? < ????; (?? ? ????)/(???? ? ????), ???? ≤ ?? ≤ ????; 1, ?? > ???? }

其中????和????分别为PCM的固相线和液相线熔化温度。在本研究中，????设定为与PCM的熔化温度????相同的值，????设定为???? + 0.1 K。??(??)为糊状函数，定义于公式(6)：

(6) ??(??) = ???? ???? [1 ? ??(??)]2 / (??3(??) + ε)

其中????为糊状系数，ε为用于避免奇异性的小值。在本研究中，???? = 1×10? kg/(m3·s)，ε = 0.001。本研究采用的????和ε的经验值也用于其他LHTES研究[13]。

边界条件定义于公式(7)。在PCM区域的所有边界上对速度场施加无滑移条件(公式7a)，对压力施加零梯度条件(公式7b)。在热源表面上，在HCM区域规定恒定的热通量q?″(公式7c)。PCM和HCM区域的所有其他外边界为绝热(公式7d, 7e)。在界面处，温度和温度梯度均假设为连续的(公式7f, 7g)。

(7a) ???? = 0, 在 Γ?,? ∪ Γ_c 上
(7b) ??·????? = 0, 在 Γ?,? ∪ Γ_c 上
(7c) ???? ??·???? = q?″, 在 Γ?,? 上
(7d) ??·???? = 0, 在 Γ?,? 上
(7e) ??·????? = 0, 在 Γ?,? 上
(7f) ???? = ???, 在 Γ_c 上
(7g) ???? ??·????? = ??? ??·????, 在 Γ_c 上

2.2 用于优化的自然对流模型

在优化过程中，根据以下假设简化了2.1节中描述的自然对流模型。

如图2(b)所示，引入设计密度场??来表示PCM区域(?? = 1)和HCM区域(?? = 0)。因此，在优化期间，PCM和HCM区域被合并为单一计算域????。为了捕捉PCM和HCM之间的材料性能变化，定义了比热、导热系数和潜热的插值函数。这些性能采用实体各向同性材料带惩罚(SIMP)方法[12]进行插值，如公式(8)所示。

(8a) C???(??) = ????????,?? + (???????,? ? ????????,??)(1 ? ??)^(????)
(8b) ??(??) = ???? + (??? ? ????)(1 ? ??)^(????)
(8c) L?(??) = ???????? ??^(????)

其中下标"p"和"h"分别代表PCM和HCM，????、????和????为惩罚因子，在本研究中均等于3。此外，在材料属性与温度无关的假设下，PCM和HCM区域中的焓可以统一表达为公式(9)：

(9) ?(??, ??) = C???(??)?? + L?(??)??(??)

为了对HCM区域中的速度进行惩罚，将Brinkman摩擦项[29]引入动量方程，定义于公式(10)：

(10) ????(??) = ???????? + (???????? ? ????????) · ????(1 ? ??)^(????) / (????(1 ? ??)^(????) + ??)

其中???????? = 0，???????? = ????/(??2????)。这里??表示特征长度，???? = 0.1为惩罚因子，????为Darcy数，本研究取为10??。

在优化过程中，评估目标函数的灵敏度需要相场??(??)的梯度。然而，公式(5)中定义的分段连续相函数在?? = ????和?? = ????处不可微。因此，在优化过程中用光滑近似替代，如公式(11)所示：

(11) ??(??) = 1 / (1 + exp(?ξlog(T ? ????)))

其中ξlog为陡度因子，本研究设为16。

因此，PCM和HCM区域的控制方程可以统一为公式(12)：

(12a) ?·?? = 0
(12b) ???/??? + ??·??? = ???? + ?·(???????) ? ????(T ? ????????)?? ? ??(??)?? + ????(??)??
(12c) ??/??? + ??·?? = ?·(??(??)???)

统一域的边界条件列于公式(13)：

(13a) ?? = ??, 在 Γ? ∪ Γ? 上
(13b) ??·??? = 0, 在 Γ? ∪ Γ? 上
(13c) ???(??)??·??? = q?″, 在 Γ? 上
(13d) ??·??? = 0, 在 Γ? 上

2.3 纯导热模型

除了自然对流模型外，在优化设计的后分析中还采用了纯导热模型。对纯导热模型和自然对流模型所得熔化行为的比较分析，提供了自然对流在熔化过程中作用的物理见解。本节详细描述纯导热模型的控制方程和边界条件。

在HCM区域，扩散方程公式(1)继续作为控制方程。在PCM区域，忽略自然对流，由于流体运动完全由自然对流驱动，因此整个域中的速度为零。因此，省略公式(3a)中的连续性方程和公式(3b)中的动量方程，公式(3c)中的能量方程简化为：

(14) ????/??? = ?·(?????????)

公式(7c)–(7g)给出的边界条件同样适用于纯导热模型中的温度场。

2.4 用于优化的导热模型

纯导热模型也用于优化过程。在3节所述的两阶段优化策略中，第一阶段使用该模型获得中间设计，该设计随后用作第二阶段（涉及自然对流模型）的初始设计。

控制方程描述如下。在与2.3节相同的区域中，消除公式(12a)、(12b)给出的连续性和动量方程，公式(12c)定义的能量方程简化为：

(15) ??/??? = ?·(??(??)???)

公式(13c)、(13d)中指定的边界条件适用于纯导热TO模型。

2.5 完整模型与优化模型的比较

在本研究中，所有数值算例（包括优化和验证算例）均使用开源有限体积软件包OpenFOAM v7.0进行。完整自然对流模型(2.1节)的CFD求解器在作者先前的研究[30], [31]中开发并经过实验验证，表明数值预测与实验测量结果吻合良好。因此，本研究将CFD模拟结果视为真实基准。

在优化中，使用带有插值材料属性且合并PCM/HCM域的简化自然对流模型(后文称为"优化NC模型")以降低计算复杂度(2.2节)。以4.1节中将详细讨论的拓扑优化设计为例，研究优化NC模型与完整NC模型之间的精度和偏差。所研究设计的几何形状如图3(a)所示，其中黑色和灰色区域分别表示HCM和PCM。在板的底部表面施加恒定热通量q?″ = 1.18 W/cm2。采用上述四个模型在该设计上进行不同保真度的模拟，图3(b–c)展示了预测的平均底板温度。完整NC模型(2.1节)和优化模型(2.2节)的预测结果分别以蓝色和红色曲线呈现在图3(b)中。模型的相对误差定义为：

(16) 相对误差 = |??full ? ??opt| / |??full ? ????|

其中??full、??opt、????分别为完整模型、优化模型预测的温度和PCM的熔化温度。优化模型与完整NC模型具有合理的吻合度，在最终时间步（即关注量）的相对误差为7.90%。类似地，完整导热模型(2.3节)和优化模型(2.4节)的预测结果如图3(c)所示。优化模型在整个熔化过程中具有较高的精度，结束时的相对误差为6.49%。

图3. 完整模型与优化模型的比较。

3. 优化方法

本节将详细讨论优化问题及其数学表述。由于非线性、局部极小值的存在以及使用自然对流模型进行优化时的高计算成本，本研究采用两阶段优化流程以获得改进的设计。两阶段策略的细节在3.1节中讨论，随后在3.2节和3.3节分别给出导热优化阶段和对流优化阶段的数学表述。

3.1 两阶段优化策略

在标准的基于密度的TO过程中，优化通常从均匀灰色设计开始，并通过连续增加Heaviside函数的陡度因子，使设计在过程中逐渐变得清晰。然而，早期设计中的显著灰度导致不利影响。具有显著灰度的设计会导致PCM和HCM区域之间边界模糊，从而导致速度预测不准确。在这些条件下，自然对流流无法被有意义的预测，因此对流换热变得不可靠。结果，表示不良的物理场可能将设计推向错误的局部极小值。此外，使用自然对流模型的前向模拟在计算上远比纯导热模型昂贵。在每个时间步，前者需要迭代求解耦合的速度-压力-温度场(约100次迭代)，而后者仅求解温度场(一次迭代)。考虑到所需的大量时间步，计算时间增加到了不实用的水平。因此，在优化早期阶段（当设计具有显著灰度时）使用自然对流模型既不具物理意义也不具计算效率。本研究采用两阶段优化策略来解决此问题。

两阶段优化过程如图1(a)所示。在第一阶段，使用2.4节讨论的纯导热模型进行TO过程，获得导热优化的设计。该过程从均匀灰色设计开始，设计逐渐演变为具有清晰边界的锐利、明确的結構。与自然对流模型相比，基于导热的优化在求解原始方程和伴随方程方面在计算上更加高效。然而，纯导热模型的较低保真度需要在第二阶段优化中对自然对流效应进行修正，从而进一步改善设计性能。

在第二阶段，使用2.2节讨论的自然对流模型进行优化。采用界面过滤结构优化(IFSO)方法[25]在整个优化过程中保持锐利边界，确保每次迭代时速度和温度场在物理上有意义。导热优化的设计将作为IFSO过程的初始设计。3.2节将证明，导热优化的设计表现优于传统设计，因此适合作为第二阶段的初始设计。通过利用初始设计的信息，第二阶段优化中的总迭代次数可以显著减少。与从灰色设计开始的常规拓扑优化通常需要300–400次迭代相比，IFSO方法仅需50–100次迭代即可实现高质量设计。迭代次数的减少使得3D几何的优化成为可能，如4.2节案例2(立方域中的3D设计)和4.3节案例3(环形域中的3D设计)所演示的。

3.2 第一阶段：纯导热拓扑优化

本节介绍第一阶段所用导热优化问题的数学表述。

对于给定持续时间????????，优化目标是最小化放电结束时热源表面上的空间平均温度，定义为：

(17) ?? = (1/|Γ?|) ∫Γ? T|????? dS

其中|Γ?|为热源表面的面积。

对原始密度场??施加具有常数半径值的Helmholtz型PDE滤波器[32]和Heaviside投影[33]，以正则化设计并消除网格依赖性问题。PDE滤波器的控制方程和边界条件列于公式(18)：

(18a) ?R2?2γ? + γ? = ??, 在 ???? ∪ ???? 中
(18b) ??·?γ? = 0, 在 Γ? ∪ Γ? 上

其中γ?为滤波后的设计密度场，通过求解上述方程获得。Heaviside投影定义为：

(19) γ? = H(γ?) = [tanh(βH·η) + tanh(βH·(γ? ? η))] / [tanh(βH·η) + tanh(βH·(1 ? η))]

其中βH和η = 0.5分别为Heaviside投影的陡度因子和阈值。在本研究的第一阶段优化中，βH初始设为1，每20次迭代增加1.25倍，直至达到32。

因此，第一阶段优化问题的数学表述可总结为：

(20a) min_?? ?? = (1/|Γ?|) ∫Γ? T|????? dS
(20c) s.t. (15), (13c), (13d), (18), (19), 0 ≤ ?? ≤ 1

在本研究的两个阶段优化过程中均未考虑HCM的体积约束。

第一阶段的整体优化工作流程如图4(a)所示。每次迭代中，原始密度场??通过公式(18)的PDE滤波器和公式(19)的Heaviside投影处理，得到投影密度场γ?。然后求解瞬态原始方程公式(15)及其边界条件公式(13c)、(13d)以获得温度场，随后评估成本函数。成本函数??对优化变量??的灵敏度????/????通过空间连续、时间离散的伴随方法进行分析，该方法已在作者先前的工作中采用。

对伴随方程详细推导感兴趣的读者可参考文献[16]和[17]。计算得到的灵敏度随后通过移动渐近线法（MMA）[34]用于更新密度场γ。优化过程以迭代方式进行，直至达到最大迭代次数。本研究中使用最大迭代次数作为停止准则，因为最大设计变化量和连续两次迭代中目标函数变化量等其他准则难以稳定实现。在第一阶段，所有算例的最大迭代次数均设为340次，以确保当Heaviside陡度因子βH每20次迭代增大1.25倍时，能够到达其最终值16。

图4. 传导拓扑优化与对流界面过滤结构优化的流程图。

3.3. 第二阶段：对流界面过滤结构优化

本节阐述第二阶段中使用的IFSO方法的数学表述。IFSO方法的整体迭代流程如图4(b)所示，其中概述了算法的主要步骤。与基于密度的拓扑优化方法类似，IFSO方法使用密度场γ来表示设计。在PDE滤波器中采用非均匀半径场R，其值作为优化变量并在每次迭代时更新。结合非均匀PDE滤波器和Heaviside投影，设计可以通过连续的边界演化进行修改造。关于PDE滤波器诱导设计边界演化机制的更多细节，可参见我们之前的工作[25]。

整个优化过程分为两个循环：用于滤波器的内循环和外循环。在内循环中，滤波器半径通过式(21)中列出的优化问题进行更新。其中R?是归一化滤波器半径，R?_min是其允许的最小值。本研究中，所有数值算例均使用R?_min = 1×10??。归一化半径R?与半径场R之间的关系由下式给出：R2 = R2_max × R?，其中R_max是允许的最大滤波器尺寸。在优化过程中，自然对流模型方程(12)和(13)控制着物理场(u, p, T)的预测，并且目标函数的灵敏度是关于优化变量R?计算的。与第一阶段相同，灵敏度使用空间连续、时间离散的伴随方法计算。关于伴随方程的推导以及关于R?的灵敏度表达式，可参见文献[25]。在正确计算灵敏度后，采用MMA方法作为优化器来更新优化变量R?。内循环迭代直至达到最大迭代次数N_R。本研究中，N_R = 5。

外循环更新设计密度场γ。每完成一轮内循环后，密度场通过以下规则进行更新：其中上标表示外循环迭代次数，下标表示内循环迭代次数。式(23)表明，前一轮外循环i的最终内步骤N_R处的投影场γ?被传递到外循环i+1中，作为原始密度场γ用于后续的内循环迭代。除了更新密度外，滤波器半径也进行了调整，以减少前一外循环的最终内循环迭代与前一外循环的初始内循环迭代之间的投影密度场之间的差异。通过求解以下优化过程获得。

根据与第一阶段相同的原因，最大迭代次数被用作对流优化阶段的停止准则。整体过程的停止准则是规定的最大迭代次数，达到该次数后优化终止。本研究中，第二阶段所有数值算例均使用最大迭代次数100次，因为这一数量足以确保目标函数近似收敛（如第4.1节所示），而不会导致过高的计算成本。

4. 数值结果

本节将第3节讨论的优化方法应用于三个数值算例。第一个算例（第4.1节）考虑一个2D矩形域，具有三种不同的高度。优化后的设计在多区域网格上进行重构，其中相变材料（PCM）区域和热传导材料（HCM）区域由不同的边界分隔。对优化后的设计以及几种常规设计进行全自然对流模型（第2.1节）的计算流体力学（CFD）模拟，以进行性能验证和比较。此外，还基于全传导模型（第2.3节）进行了模拟，以与全自然对流模型结果进行比较，并说明自然对流对熔化行为的影响。在第二个算例（第4.2节）中，优化在3D立方体域上进行，而在第三个算例（第4.3节）中，考虑一个水平放置的3D环形域。与2D算例类似，通过全自然对流CFD模型对优化后的3D设计进行数值验证，并将其热性能与传导优化设计进行比较。在两种情况下均观察到显著改善。

优化和CFD模拟基于开源有限体积（FV）软件OpenFOAM v7.0实现。在所有数值算例中，冰水用作PCM，铝锌合金7075-T6用作HCM。相关材料特性见表2。本研究中使用的水的热膨胀系数对应于T = 25°C时的水。如下文所示，本文所考虑算例中，最终时刻的温度大致在20-25°C范围内。由于优化目标集中在最终时间步长的温度上，因此选择相应温度范围内的材料特性，使优化过程能够更好地针对期望的热行为。时间离散化采用一阶欧拉格式进行。空间上，对流项采用二阶迎风格式，而扩散项采用二阶线性格式。速度-压力耦合通过压力耦合方程的半隐式方法（SIMPLE）处理。优化算例在一台双路AMD EPYC 7K62（96核CPU）工作站上运行，3D算例在12个处理器上并行化。传导优化阶段（第一阶段）的计算时间约为12小时，而对流优化阶段（第二阶段）约为3-7天。

表2. 所用PCM和HCM的材料特性

| 材料特性 | 水（PCM） | Al 7075-T6 |
| :--- | :--- | :--- |
| 密度 [kg/m3] | 1000 | 2613 |
| 比热容 [J/kg·K] | 4180 | 883 |
| 导热系数 [W/m·K] | 0.613 | 20 |
| 潜热 [J/kg] | 334 000 | – |
| 熔点 [°C] | 0 | – |
| 黏度 [Pa·s] | 0.001 | – |
| 膨胀系数 [1/K] | 0.000257 | – |

4.1. 算例1：2D设计

本算例中，几何体由一个尺寸为40 mm长、20 mm高的2D矩形域（设计域）和位于设计域正下方的16 mm × 1 mm底板（非设计域）组成，如图5(a)第一张图所示。底板底面施加了均匀热流密度q'' = 1.18 W/cm2，对于76.2 mm的非面内深度，对应加热器功率为14.4 W。初始温度设为-1°C。放电持续时间的确定如下。使用以下方程计算完全熔化所需的理想时间：t_ideal = ρ_p × L_p × V_Ωd / Q?，其中V_Ωd是设计域的体积，Q?是放电功率。根据式(25)计算的理想时间基于整个域内均匀升温以及所有输入能量转化为潜热的假设。在本算例中，理想时间计算为t_end = 1414 s。在实际熔化过程中，温度在域内并非均匀升高，且只有部分输入热量以潜热形式储存，其余部分以显热形式储存在PCM和HCM中。因此，实际熔化时间t_end长于该理想估计值。本研究中，t_end随机选取在理想熔化时间t_ideal的1.0-1.15倍范围内。据此，本算例的实际结束时间设为t_end = 1600 s（1.1 × t_ideal）。利用对称性，仅研究域的一半以减少计算成本。计算域通过100 × 100的2D结构化FV网格离散化，时间域通过时间步长Δt = 2 s离散化。

在传导优化阶段（第一阶段），PDE滤波器半径设为R = 1.5个单元。如第3.2节所述，Heaviside阈值设为η = 0.5。Heaviside投影的陡度因子初始设为β_H = 1，并每20次迭代增大1.25倍，直至达到β_H = 32。第一阶段总迭代次数为340次。传导拓扑优化过程的设计演化如图5(a)所示。优化从γ = 0.7的均匀灰色设计开始。随着优化进行，翅片结构出现，灰色区域逐渐减少。优化结束时获得清晰的设计。在优化后的设计（称为传导拓扑优化设计）中，在左上和右上角形成了分支翅片。此外，在左下和右下角，水平翅片从主分支翅片伸出。图5(b)显示了目标函数值的演化，其在优化过程中平滑收敛。

在进行第二阶段优化之前，首先通过全自然对流CFD模型对传导拓扑优化设计进行数值模拟，并与4种常规设计进行比较。在CFD模拟中，时间步长减小到Δt = 0.1 s以提高精度。本研究中，考虑了三种不同翅片高度的平行翅片（PF）设计和一种分形翅片设计，如图6(a)所示。所有四种设计的HCM体积与传导拓扑优化设计相同，以确保公平比较。CFD模拟预测的源面平均温度如图6(b)所示。在设计条件下（t = 1600 s），六翅片PF设计（PF-6）在四种常规设计中表现出最低温度，平均温度为22.28°C。与PF-6相比，传导拓扑优化设计进一步将温度性能改善了0.61°C（2.4%）。结果表明，即使在优化过程中不考虑全自然对流模型，传导拓扑优化设计在包含自然对流的后CFD验证中也优于常规设计。因此，传导拓扑优化设计被认为是第二阶段优化的合适初始设计。PF-6、分形和传导拓扑优化设计的物理场如图6(c)所示，其中域的左半部分表示相场，右半部分显示速度场。在PF-6设计中，相变优先发生在下部区域。在1200 s时，底部区域完全熔化，而大量上部区域的PCM仍保持固态。相比之下，分形翅片设计促进了上部区域更快的熔化。在1400 s时，顶部区域完全处于液相，而左下角仍保持固态。在传导拓扑优化设计中，与另外两种设计相比，上下区域的熔化速度更加均衡。在1400 s时，传导拓扑优化设计在顶部区域的未熔化PCM少于PF-6设计，在底部区域的未熔化PCM少于分形设计。更均匀的熔化行为导致传导拓扑优化设计具有更低的底板温度。

第二阶段优化过程如图5(c-d)所示。在IFSO算法中，最大PDE滤波器尺寸R_max设为1.5个单元，与第一阶段一致。Heaviside函数的陡度参数β_H在整个优化过程中保持为16不变，以维持尖锐边界。

设计密度场??每5个内循环（????=5）更新一次，总迭代次数为100次。图5(c)表明，在第二阶段优化过程中，设计始终保持清晰的边界。顶部分支翅片逐渐减少，底部侧翅片逐渐缩短。

下载：下载高分辨率图片（846KB）下载：下载全尺寸图片

图7. 40 mm × 20 mm案例中传导TO设计与对流传热IFSO设计的CFD模拟结果。

本案例中第一阶段的计算时间列于表3的第一行，第二阶段列于第二行。每次优化迭代中，传导阶段仅需1.22分钟，而自然对流阶段则需要39.26分钟。这一显著差异表明，自然对流模型的计算成本远高于传导模型。

表3. 采用传导模型和自然对流模型的优化过程计算时间。

| 案例 | CPU数量 | 总CPU时间 | 每次迭代CPU时间 |
|------|---------|-----------|----------------|
| 2D矩形，传导 | 1 | 416分钟 | 1.22分钟 |
| 2D矩形，对流 | 1 | 3926分钟 | 39.26分钟 |
| 3D矩形，传导 | 12 | 523分钟 | 1.54分钟 |
| 3D矩形，对流 | 12 | 5193分钟 | 51.93分钟 |

第二阶段中目标函数的收敛历史如图5(d)所示。目标函数总体呈下降趋势，但观察到的跳跃是由于频繁更新密度所致。第二阶段优化后的设计（记为对流传热IFSO）同样通过完整的自然对流CFD模型进行了数值验证，其性能如图7(a)所示。将自然对流纳入优化过程后，对流传热IFSO设计在设计时刻的平均温度相比传导TO设计进一步降低了2.54°C。温度性能的百分比改善可由公式(26)计算：

$$\text{改善} = \frac{T_{cond|tend} - T_{conv|tend}}{T_{cond|tend} - T_m} \times 100\%$$

其中$T_{cond}$和$T_{conv}$分别代表传导优化TO设计和对流优化IFSO设计的平均基底温度。根据该公式，本案例中的改善率为11.7%。与最佳传统设计相比，性能改善了3.15°C（14.1%）。传导TO设计和对流传热IFSO设计的相场和速度场如图7(b)所示。尽管上方区域没有分支翅片，对流传热IFSO设计在该区域仍能保持相当的熔化速率。例如，在1200 s和1400 s时，上方区域的熔化前沿与传导TO设计中观察到的相似。通过将那些HCM翅片替换为PCM，对流传热IFSO设计可以储存额外的潜热能量，从而降低放电过程结束时的温度。

为了阐明自然对流的物理影响，传导TO设计也使用完整的传导模型（见2.3节）进行了模拟，以便进行对比。温度场和熔化前沿与自然对流CFD模型的预测结果进行了比较，如图8所示。在图8(a,b)中，域的左半部分代表相场，右半部分显示温度场。第一行对应仅传导预测，第二行代表对流预测。在第三行中，左侧图形显示了从相场中提取的熔化前沿，作为固体和液态PCM之间的界面，蓝色和红色线分别表示仅传导和对流预测。右侧面板展示了来自对流模型的速度场。在早期阶段（300 s），传热由传导主导，自然对流的影响可以忽略不计，正如低速度幅值所表明的。在此阶段，两个模型预测的熔化前沿几乎相同。当形成足够厚的液态PCM层且温度梯度变大时，自然对流开始影响熔化行为。自然对流驱动热液态PCM向上运动，通过对流运动将热量从底部输送到顶部区域。因此，顶部区域的熔化速率加快，而底部区域则减慢。在900 s时，可以观察到自然对流模型预测的熔化前沿（红线）在中间翅片附近开始领先于传导模型预测的熔化前沿（蓝线）。这种"向上驱动"机制在后期变得更加明显。在1440 s和1500 s时，自然对流模型中顶部区域的PCM完全熔化，而传导模型仍然预测存在固态PCM。相反，自然对流模型预测的熔化前沿相比传导模型有所延迟。这清楚地表明，自然对流加速了顶部区域的熔化速率，同时减缓了底部区域的熔化速率。因此，在考虑自然对流的设计中，应减少或消除顶部的HCM结构，以平衡顶部和底部区域之间的熔化速率，从而促进设备内更均匀的热量分布和储存。

下载：下载高分辨率图片（680KB）下载：下载全尺寸图片

图8. 仅传导和对流CFD预测的温度、相场和速度场。（对于图例中颜色引用的解释，读者请参阅本文的网络版本。）

除了上述20 mm高度案例外，还在相同优化设置下研究了30 mm和40 mm高度的两个附加案例。由于对称性，同样仅考虑域的一半。计算域分别采用100 × 150网格（30 mm案例）和100 × 200网格（40 mm案例）进行离散化。两个案例的优化结果分别如图9(a,b)所示。两个案例中均观察到类似的形状演化趋势。在传导TO阶段，优化器倾向于在整个设计域内分布HCM，在远场形成高耸的分叉翅片。在对流IFSO阶段，这些高耸的分叉翅片逐渐被移除，以平衡上下区域之间的熔化速率。在30 mm案例中，阻碍核心区域对流的中央翅片被进一步消除，允许底部中心的热流体上升至顶部区域。两个对流传热IFSO设计均通过完整的自然对流模型进行了验证，其对应的传导TO设计也一并进行了验证。CFD模拟结果如图9(c,d)所示。两个案例均观察到性能改善。对于30 mm高度案例，设计时刻（2400 s）的平均温度降低了2.43°C（10.6%）；对于40 mm高度案例，降低了2.93°C（15.0%）。在所有三个案例中，自然对流都被发现增加了顶部区域的熔化速率，同时降低了底部区域的熔化速率。对流传热优化器相应地移除了上部翅片以平衡熔化速率，从而提高了设计时刻的温度性能。

下载：下载高分辨率图片（588KB）下载：下载全尺寸图片

图9. 2D 40 mm × 30 mm和40 mm × 40 mm案例的设计优化和CFD验证。

下载：下载高分辨率图片（1MB）下载：下载全尺寸图片

图10. 3D矩形案例的设计优化和CFD验证。（对于图例中颜色引用的解释，读者请参阅本文的网络版本。）

**4.2. 案例2：立方域中的3D设计**

本研究采用尺寸为40 mm × 40 mm × 20 mm的3D矩形域作为设计空间，其中20 mm维度沿重力方向对齐。设计由底部一块16 mm × 16 mm × 1 mm的底板支撑。在底板上施加恒定7.6 W的加热器（相当于$\dot{q}'' = 2.97 \text{ W/cm}^2$）。根据公式(25)，理想熔化时间为$t_{ideal} = 1406$ s，实际总放电时间指定为理想熔化时间的约1.05倍，对应$t_{end} = 1520$ s。初始温度设为-1°C。利用水平和深度方向的对称性，域被分为四个相同的部分，仅考虑其中四分之一作为计算域。采用50 × 50 × 50的结构化网格对计算域进行离散化，时间步长为$\Delta t = 2$ s用于时间离散化。在传导优化阶段（第一阶段），采用与4.1节相同的PDE滤波器半径（$R = 1.5$ cell）和Heaviside投影方案。参数$\beta_H$初始设为1，每20次迭代增加1.25倍，直至达到$\beta_H = 32$。第一阶段总迭代次数为340次。在对流IFSO阶段（第二阶段），PDE滤波器的最大半径$R_{max}$也设为1.5 cell，Heaviside陡峭度固定为$\beta_H = 16$。设计密度场每5次迭代更新一次，第二阶段总迭代次数为100次。

整个设计演化历史如图10(a)所示。在传导优化阶段，底板上生成了九个根部翅片，排列成大致3 × 3的模式。在根部翅片上，顶部和两侧形成了分叉翅片，以促进热量向远离中央热源区域的PCM传导。值得注意的是，传导TO设计的结构与文献[24]中报告的优化3D设计类似。在对流IFSO阶段（第二阶段），上部翅片结构——特别是复杂的分叉翅片——在优化过程中逐渐被移除，而下部翅片结构则大体保留。为了更好地可视化和比较，优化后的传导TO和对流IFSO设计分别如图10(b,c)所示。从水平视图（图10(b,c)第一行）可以清楚地观察到，IFSO设计的顶部结构被减少，导致整体结构高度低于传导TO设计。

传导阶段和自然对流阶段的计算时间分别列于表3的第三行和第四行。传导模型每次迭代需要1.54分钟，而自然对流模型需要51.93分钟。这一巨大差异凸显了自然对流阶段的计算成本之高，这也证明了所提出的两阶段方法的必要性。

与4.1节类似，基于完整的自然对流模型进行CFD模拟以验证传导TO和对流IFSO设计。在CFD模拟中，时间步长为0.1 s。性能如图10(d)所示。在早期熔化阶段，对流传热IFSO设计的平均基底温度高于传导TO设计，因为早期传热由传导主导，而对流IFSO设计含有较少的HCM结构用于热传导。在后期阶段，一旦形成足够厚的液态PCM层，对流成为主导的传热机制，对流传热IFSO设计的表现优于传导TO设计。在设计时刻（1520 s），基底温度降低了2.61°C，代表13.2%的改善。

自然对流的强度可以通过PCM区域的最大速度来表征[19]。在本研究中，两个设计的最大速度如图10(e)所示。在大部分放电期间，对流传热IFSO设计中的速度高于传导TO设计。在峰值速度时刻，对流传热IFSO设计中的速度幅值约为传导TO设计的2.1倍，表明对流传热IFSO设计中自然对流强度增强。这可由对流传热IFSO设计中更高的温度梯度来解释。从公式(3b)中，其中自然对流项由$\beta_f(T_p - T_{ref})\mathbf{g}$给出，可以看出流速是由温度梯度驱动的。因此，对流传热IFSO设计中较大的温度梯度导致更高的对流强度。

还对两个设计的熔化行为进行了研究。传导TO和对流IFSO设计在全域四分之一范围内的熔化前沿分别如图10(f,g)所示。蓝色区域代表固态PCM。

在 400 秒和 800 秒时，导热 IFSO 设计的顶面（由黄色和绿色方块标记）出现孔洞，表明这些区域的相变材料（PCM）已经熔化。与此同时，对流 IFSO 设计顶面仍保留着大片固态 PCM 区域，表明对流 IFSO 设计中顶部的熔化速度比导热 TO 设计慢。对流 IFSO 设计的顶部区域在 1200 秒时才开始熔化（由红色方块标记）；然而，在同一时间，导热 TO 设计的该区域已几乎完全熔化。为了量化熔化速度，将域沿重力方向划分为两个相同的子域。计算了每个子域中的液态分数，并在图 10(h, i)中显示，分别对应上半域和下半域。从图 10(h) 可以明显看出，对流 IFSO 设计的一半在整个过程中表现出降低的熔化速度，导致完全熔化时间大幅延迟。作为补偿，下半域经历了加速的熔化率，如图 10(i) 所示。

下载：下载高分辨率图像 (1MB) 下载：下载全尺寸图像 图 11.3D 环形案例的设计优化与 CFD 验证。（关于图注中颜色引用的解释，请参阅本文的网络版本。）

4.3. 案例 3：环形域中的 3D 设计

本示例考虑了一个重力垂直作用的水平环形区域。如图 11(a) 所示，环形区域的内半径和外半径分别为 3 mm 和 25 mm。在域内，从 ??=3mm 到 5 mm 的区域是基板，填充有高导热材料（HCM），而外域 ??=5?25mm 是设计域。环形区域的非平面深度为 40 mm。在环形区域的内表面施加了一个总功率为 15.7 W（对应 2.09 W/cm2）的加热器。根据公式 (25) 计算的理想熔化时间为 1600 秒。在本示例中，实际放电时间规定为理想熔化时间的 1.0 倍，对应 ??????????=1600s。初始温度仍为 ?1°C。沿图 11(a) 中的虚线将域二等分，仅计算一半以减少计算成本。空间域使用结构化网格进行离散化，径向、方位向和展向的控制体数量分别为 76、60 和 40 个。优化过程中的时间步长为 2 秒。两个优化阶段的参数与 3D 矩形案例（第 4.2 节）中的参数保持相同，除了第一阶段中的 PDE 过滤尺寸 ?? 和第二阶段中的最大过滤半径 ???????? 被设置为 2 个单元。

优化过程如图 11(b) 所示。在导热优化阶段，逐渐生成翅片以促进热量从内表面传导至外表面。这些翅片源自基板，并在靠近外表面处分叉。在第二阶段，将对流 IFSO 方法应用于优化的导热 TO 设计。可以观察到，在 IFSO 过程中，上半部分的翅片结构被简化，部分最终被移除，而下半部分的翅片结构则大部分得以保留。两个阶段的优化设计分别如图 11(c, d) 所示。两个设计的正视图更清晰地展示了它们的几何差异。在对流 IFSO 设计中，上区域的翅片比导热 TO 设计中的更粗糙且更简单，靠近外表面的分叉翅片也减少了。在底部区域，两种设计中的翅片结构相似。这种形状演变趋势与第 4.2 节中所示的矩形案例一致。

使用完整的自然对流模型对两种优化设计进行了 CFD 验证。CFD 模拟预测的热源（内表面）平均温度绘制在图 11(e) 中。在设计时间（1600s）时，导热 TO 设计的温度为 16.92°C，而对流 IFSO 设计的温度为 15.24°C。对流 IFSO 设计实现了 1.68°C 的温度降低，相当于 9.9% 的性能提升。PCM 场中的最大速度值如图 11(f) 所示。在大多数放电过程中，对流 IFSO 设计中的速度高于导热 TO 设计中的速度。在峰值速度时，对流 IFSO 设计的速度大小是对流 TO 设计的 1.29 倍。

两种设计的熔化过程如图 11(g, h) 所示，其中蓝色区域表示未熔化的 PCM。在 600 秒和 900 秒时，导热 TO 设计的外表面出现孔洞，而对流 IFSO 设计外表面附近的 PCM 仍保持固态，表明在对流 IFSO 设计中，上区域的熔化过程被减缓。在后期阶段，随着自然对流成为主导因素，对流 IFSO 设计中的熔化率开始加速。在结束时间（1600 秒）时，可以观察到两种设计中的上区域均已完全熔化。与矩形案例类似，计算域通过水平面分为上半部分和下半部分。两半部分中定量的潜热储存如图 11(i, j) 所示。在域的上半部分，对流 IFSO 设计中的潜热储存率减速，导致完全熔化时间延迟，而在下半部分，作为补偿，早期熔化率得到增强。这些结果表明，对流 IFSO 阶段中的形状演变缓解了由浮升力驱动的上向流动引起的上下熔化速率差异，从而增强了热性能。

5. 结论

本研究提出了一种用于潜热储热（LHTES）设备拓扑优化（TO）的两阶段优化方法，解决了一个涉及自然对流和相变的瞬态热 - 流问题。在第一阶段，使用仅导热模型进行优化。尽管第一阶段忽略了自然对流等关键物理现象，但导热优化设计仍表现出强劲的性能，优于传统设计。同时，计算成本大幅降低（导热模型中每时间步 1 次迭代 vs. 对流模型中每时间步≈100 次迭代）。在第二阶段，将界面过滤结构优化（IFSO）方法应用于第一阶段获得的 TO 设计。采用包含自然对流的完整模型来细化设计并进一步提升其性能。借助第一阶段生成的良好初始设计，第二阶段优化迭代次数可从 300–400 次减少到 50–100 次。两阶段方法加速的优化程序和降低的计算成本使得 3D LHTES 设备的拓扑优化成为可能，否则这在计算上是不可行的。

对 2D 和 3D 示例进行了两阶段优化，并通过计算流体力学（CFD）模拟对优化设计进行了后验证。在 2D 示例中，第一阶段导热优化设计四种传统设计进行了比较，包括具有不同翅片数量的平行翅片配置和分形翅片设计。导热 TO 设计比最佳传统设计高出 2.4%。此外，第二阶段的对流 IFSO 设计相比导热 TO 设计性能提升高达 15.0%。还研究了 3D 立方体域和环形域，性能提升高达 13.2%。观察到一种“驱动向上”机制，并用于解释自然对流在设计优化中的作用，即浮升力驱动的流动将热量从底部带至顶部，增强上区域的熔化，同时减缓下区域的熔化。因此，在优化设计中应减少上部高导热结构，以补偿不均匀的熔化行为。在所有 2D 和 3D 示例中均观察到相应的形状变化。

值得强调的是，本文提出的两阶段策略并不局限于 LHTES 应用。它构成了多物理场耦合问题的一种通用范式，其中第一阶段采用简化的物理模型，而第二阶段调用完整物理模型。这样做可以显著降低计算开销和运行时间，使大规模优化成为可能。

总之，本研究的关键观察结果可总结如下：
•导热优化设计优于优化后的传统设计，从而为基于对流的 IFSO 阶段提供了有效的初始设计。
•通过考虑自然对流，可进一步提升热性能。与导热优化设计相比，基于对流的 IFSO 设计在温度性能上实现了高达 15.0% 的提升。
•自然对流产生浮升力驱动的流动，通过对流将热量从底部传输到顶部，增强上区域的熔化，同时延缓下区域的熔化。因此，经过调整以缓解由此产生的熔化速率失衡的设计实现了热性能的提升。
•对于实践者而言，本研究的观察结果表明，当自然对流效应显著时，对称结构并不是最优的。减少顶部翅片结构并增强底部翅片结构是提高热性能的有用策略。

CRediT 作者贡献声明

Tianye Wang：撰写 – 初稿，验证，软件，方法论，概念化。
Allison Mahvi：撰写 – 审阅与编辑，监督，概念化。
Xiaoping Qian：撰写 – 审阅与编辑，监督，概念化。

Tianye Wang | Allison Mahvi | Xiaoping Qian
美国威斯康星大学麦迪逊分校机械工程系，美国