时域有限元方法与超材料中任意形状热聚集体的数值模拟

《Computers & Mathematics with Applications》：Time domain finite element method and numerical simulation of arbitrarily shaped thermal concentrators in metamaterial

　　摘要
为了实现热能的高效收集和利用，本文基于变换热力学原理推导了任意形状热集中器的材料参数。随后，为超材料区域内的控制方程建立了时域有限元公式，并对数值算法进行了全面的稳定性分析和误差估计。最后，通过数值模拟实验验证了推导出的材料参数的正确性以及所提出数值方案的有效性。

引言
早期关于超材料的研究主要集中在电磁超材料上，特别是左手材料的研究是这一领域的重要发展。所谓“左手材料”，是因为当介电常数和磁导率都为负时，电场、磁场和波矢满足左手规则。早在1968年，V. G. Veselago就基于麦克斯韦方程和材料本构关系，首次对左手材料及其异常特性进行了理论分析和预测[1]。然而，由于自然界中不存在左手材料，这一领域长期没有实质性进展或突破。直到20世纪末，Pendry等人提出周期性排列的金属杆和裂环谐振器可以分别产生负介电常数和负磁导率[2]、[3]。2000年，D. R. Smith等人基于Pendry的理论基础，通过有序排列金属杆和裂环谐振器，制造出了世界上第一个一维左手材料[4]。2002年，基于相同的设计概念，开发出了第一个二维左手材料，并通过实验确认了其负折射特性[5]。自那时起，左手材料迅速成为电磁学乃至更广泛科学界的高度活跃的研究焦点。许多看似违背经典电磁理论且与传统理解不符的奇异电磁特性得到了验证，例如完美透镜效应[6]和 backward-wave 特性[7]。在治疗电磁超材料的研究中，使用坐标变换方法来设计任意控制电磁场分布的超材料电磁器件已成为全球范围内的一个重要研究方向。这一方法起源于Pendry等人2006年在《Science》杂志上发表的论文，该论文介绍了一种通过坐标变换设计的圆柱形隐身装置[8]。此外，由于其他基本物理方程与麦克斯韦方程之间的相似性，坐标变换方法迅速扩展到了电磁学[9]、[10]、[11]、声学[12]、[13]、弹性波[14]、[15]、物质波[16]、[17]以及热力学[18]、[19]、[20]、[21]、[22]、[23]等领域。

2008年，Fan等人[24]首次将变换光学扩展到热力学领域，建立了稳态条件下的变换热力学理论，并提出了热超材料的概念，从而实现了对热流的操控。2012年，Guenneau等人[19]证明了二维瞬态热传导方程的形式不变性，推导出超材料区域内的材料参数表达式，进一步发展了变换热力学，并展示了圆柱形隐身器和集中器的设计过程和性能。此后，热超材料的研究持续进步，提出了多种能够操控热流的超材料器件。在众多热能器件应用中，热集中器作为一种有效的热能收集和捕获手段[25]，在太阳能收集、高灵敏度温度检测和热伪装等领域展现出巨大潜力。与热隐身器（将热流重新导向物体周围）不同，热集中器旨在通过调整热导率来将热流集中到特定区域，同时保持外部温度分布不变，这为提高热能收集效率提供了有希望的应用方案。2011年，Yu等人[26]首次使用坐标变换理论设计了稳态热集中器。随后，在2012年，Guenneau等人[19]将热集中器扩展到瞬态范畴，并通过模拟成功验证了其有效性。2014年，Li等人[27]基于变换热力学推导出了任意形状热集中器的热导率通用表达式。2015年，Lan等人[28]首次提出了同时具备热流集中和电流集中功能的器件，该器件由各向同性材料（铝合金和环氧树脂）组成，并通过扇形结构实现各向异性。2016年，Shen等人[29]在扇形复合多层结构中引入了形状记忆合金，通过形状记忆合金的选择性温度依赖性效应实现了热隐身和热集中功能的自动切换。2018年，Liu等人[30]提出了一种可调热超材料器件，该器件由热变换单元和旋转方形单元组成，能够在热隐身和热集中功能之间切换。然而，据我们所知，关于任意形状瞬态热集中器的文献仍然很少。此外，由于大多数先前研究都集中在稳态条件下的规则对称形状（如圆柱形）上，这些方法不适用于设计边界无法用数学函数表达的不规则形状的热集中器。

作为解决偏微分方程的最有效数值方法之一，有限元方法在处理复杂几何形状和各向异性材料特性方面具有显著优势[31]、[32]、[33]、[34]，使其成为模拟超材料器件性能的首选工具[35]、[36]、[37]、[38]、[39]。目前关于热集中器的研究主要集中在稳态条件和规则[40]、对称边界轮廓上。为了克服这些限制，本研究开发了用于分析和模拟任意形状热集中器的时域有限元方法。因此，在本文中，我们推导了两种类型任意形状热集中器的材料参数的数学表达式，并展示了瞬态条件下控制方程的有限元数学分析和数值结果。据我们所知，这是首次使用有限元方法分析和解决任意形状热集中器的工作。

本文的其余部分结构如下：第2节推导了考虑已知和未知边界轮廓条件的任意形状热集中器的材料参数。第3节对超材料区域内的时域控制方程进行了有限元稳定性估计和误差分析。第4节通过数值实验验证了推导出的材料参数的正确性和所提出数值方案的有效性，并阐明了热流从瞬态行为到稳态的物理过程。最后，第5节总结了全文。

模型控制方程和材料参数推导
众所周知，在可用的自然能源中，热能是最丰富的能源之一。由煤炭、石油和天然气等各种能源产生的热能可以保障我们的生活质量。因此，热能也成为最重要的能量耗散方式。另一方面，随着新兴技术的迅速发展和应用，器械的功耗持续增加，热能的产生也随之增加。

数值方案及其分析
在这里，我们考虑具有以下各向异性材料参数的抛物线方程：
$$
\chi_c \mathbf{T}_t - \mu \cdot (\mathbf{k}_5 (\mathbf{x}) \mathbf{\mu} \mathbf{T} = \phi(\mathbf{x}, \mathbf{t}), \quad (\mathbf{x}, \mathbf{t}) \in \Omega \times (0, T_1), \quad \mathbf{T}_2(\mathbf{x}, 0) = \mathbf{T}_0(\mathbf{x}), \quad \mathbf{x} \in \Omega, \quad \mathbf{T}_2(\mathbf{x}, \mathbf{t}) = \mathbf{g}, \quad (\mathbf{x}, \mathbf{t}) \in \Gamma_{\mathbf{D}} \times (0, T_1), \quad \nabla \mathbf{k}_5(\mathbf{x}, \mathbf{T}_2(\mathbf{x}, \mathbf{t}) \cdot \nabla \mathbf{k}_\nu = \mathbf{q}, \quad (\mathbf{x}, \mathbf{t}) \in \Gamma_{\mathbf{N}} \times (0, T_1),
$$
其中 $\Gamma_{\mathbf{D}} \cup \Gamma_{\mathbf{N}} = \Gamma = \nabla \Omega$。在后续模拟中，我们假设域的上下边界是绝热的（即 $\mathbf{q} = 0$），并且 $\Omega$ 的左右边界具有固定温度。这里，我们假设 $k(x)$ 是一个 $2 \times 2$ 的对称正定矩阵，对于任何给定的 $x \in \Omega$，有 $k^T x \geq 0$。设 $\lambda_1$ 和 $\lambda_2$ 是...

数值测试
在本节中，我们提供了一些数值结果来证明我们的理论分析的正确性和任意形状热集中器的热集中性能。我们采用了方程（49）和（50）中的完全离散有限元方案，并使用线性拉格朗日元素。

示例 1
在这个例子中，我们测试了半离散方案（33）和（34）以及完全离散方案（49）和（50）的最优收敛阶数。我们选择的精确解是 $u(\mathbf{x}, \mathbf{y}, \mathbf{t}) = e^{-t \sin(3\pi x)}$。

结论
在本文中，我们首先推导出了任意形状热集中器的理想材料参数。然后，我们开发了超材料区域控制方程的时域有限元方法，并对所提出的半离散和完全离散方案进行了稳定性分析和误差估计。数值实验验证了材料参数的准确性和数值方法的有效性。在未来的研究中，我们计划进一步探讨...

利益冲突声明
作者声明没有利益冲突。

Bin He | Haohao Han | Jiangang Zhang | Meng Chen
甘肃复杂系统分析与控制基础研究中心，兰州交通大学数学与物理学院，中国甘肃省兰州市，730070