利用高精度Galerkin无网格方法、移动克里金插值和kd树对各向异性功能梯度材料（FGM）进行瞬态热分析

《COMPOSITE STRUCTURES》：Transient thermal analysis in anisotropic FGMs by high precision Galerkin meshless method with moving Kriging interpolation and kd-tree

　　摘要：功能梯度材料（FGMs）的成分在空间中连续变化，它们在高温环境中有着重要的应用，包括航天器的热防护、核反应堆和燃气轮机叶片。然而，这些材料的不均匀性质给热分析的精确数值模拟带来了巨大的挑战。本文提出了一种高精度的Galerkin无网格方法，该方法结合了移动Kriging插值和kd树（k维树）技术，用于解决各向异性FGMs的瞬态热分析问题。所提出的方法通过移动Kriging插值构建形状函数，并明确满足了Kronecker delta性质。为了加速自适应邻域查询，引入了kd树来进行节点搜索，从而实现快速搜索操作。通过Galerkin弱形式建立了数值离散公式，并采用隐式方案进行时间离散化。通过六个例子验证了该方法的有效性，包括将数值结果与解析解以及其他数值方法的结果进行比较。验证结果表明，该方法在各向异性FGMs的瞬态热分析中具有良好的准确性、收敛性和稳定性。

引言：在航空航天推进系统、核能设备和高功率电子设备等前沿领域，材料通常需要同时承受极端的热梯度、复杂的机械载荷和恶劣的环境侵蚀[1]。在这些应用中，主要的结构部件由薄板、穿孔板和三维固体组成，它们都受到瞬态热载荷的影响，而准确预测这些载荷对于安全和可靠的设计至关重要。以可重复使用的航天器为例，其热防护系统在大气再入过程中要承受超过2000°C的空气动力加热，同时需要保持轻量级的结构设计和长寿命，这对材料和机械性能提出了几乎矛盾的要求[2][3]。功能梯度材料（FGMs）通过其在空间方向上的连续成分转换独特的设计理念，为这些挑战提供了一种突破性的解决方案[1]。这种材料架构在热障涂层、耐磨部件和生物医学植入物中展现了不可替代的优势，并被《自然材料》杂志评为未来十年最具变革性的先进材料范式之一[4][5]。在这些热和机械服务环境中，还广泛研究了功能梯度板和壳体的机械行为，包括在热和压缩载荷共同作用下的屈曲、自由弯曲振动以及更高阶的壳体理论[6][7][8]。自从Koizumi首次介绍FGMs的研究以来，复合材料领域的文献中出现了大量关于功能梯度板和壳体的综述，这些综述全面调查了热建模技术和解决方法[10]，这表明了对这类结构问题的持续研究兴趣。准确预测和优化这类异质材料在瞬态热载荷下的响应行为至关重要。传统的均匀材料分析理论和设计标准难以直接应用，因为FGMs的参数随空间连续变化，导致控制方程缺乏解析解。FGMs的性能对局部梯度非常敏感，这种敏感性源于其在服务环境中的性能定制和多场耦合，这对当前的研究构成了重大挑战。因此，迫切需要能够解决具有复杂几何形状的FGM板和三维固体中瞬态温度场的数值方法。有限元方法（FEM）一直是结构热分析的主要数值工具。对于FGMs，最常见的策略是根据局部成分在每个元素上分配有效的材料属性，使用等参分级元素[11]。更先进的FEM公式包括增强的有限元素，可以在粗网格上捕捉剧烈的热梯度[12]，以及不兼容的分级有限元素，可以在不增加自由度的情况下提高精度，这一点已在FGM固体的热传递和热应力分析中得到验证[13]。同时，基于网格的方法有一个几何要求：网格既要符合域边界，又要适应材料属性的空间变化。当考虑具有不规则边界的三维FGM体时，生成高质量的符合网格非常耗费资源，而在剧烈梯度下的网格畸变或自适应细化可能会降低解的精度，需要重新网格划分。

近年来，无网格方法已被用于异质热分析。Belytschko等人提出的无元素Galerkin（EFG）方法通过移动最小二乘（MLS）近似构建形状函数[14]。Liu等人开发了再生核粒子方法，通过修正的核函数恢复完整性[15]。Atluri和Zhu引入了无网格局部Petrov-Galerkin（MLPG）方法，通过在小的子域上使用局部弱形式来减少总体计算负担[16]。Liew等人提供了关于层压板和功能梯度板的无网格方法的全面综述[17]。在FGMs的瞬态热分析特定背景下，几种无网格方法取得了显著的成功。Sladek及其同事结合拉普拉斯变换技术，使用MLPG方法解决了三维各向异性FGM固体的瞬态热传导问题[18][19]。Zhang等人将EFG方法应用于各向异性材料的瞬态热传递[20]。Khosravifard等人开发了一种改进的径向点插值方法，用于包含内部热源的FGMs的非线性瞬态热传导[21]。在平滑粒子流体动力学（SPH）框架下，最近提出了一个用于功能梯度结构热弹性分析的对称公式[22]。由于EFG和MLPG中的MLS形状函数缺乏Kronecker delta性质，因此必须通过拉格朗日乘数或惩罚函数间接施加关键边界条件。尽管这些技术有效且应用广泛，但它们引入了额外的自由度[23][24][25]。许多无网格方案依赖于对所有节点的邻居搜索，对于涉及数万个节点的大型三维问题来说，这变得非常耗时[26]。本工作的核心贡献是将Gu的移动Kriging插值（MKI）[27]与Bentley的kd树索引[28]结合，并采用统一的Galerkin弱形式。MKI构建的形状函数本质上满足Kronecker delta性质，因此可以直接施加关键边界条件，而无需额外的乘数或惩罚；kd树通过递归二分空间划分降低了预期的邻居查询复杂度。

以下各部分的安排如下：第2节介绍了二维和三维各向异性FGMs瞬态热分析的控制公式。第3节围绕形状函数构建、基于kd树的节点自适应搜索以及时间和空间离散化推导出Galerkin弱形式。第4节通过六个数值例子验证了所提出方法的准确性、稳定性和效率。第5节给出结论。

瞬态热分析问题的控制公式：考虑一个开放有界域Ω?R^d=2,3。在这个域中，具有空间变化导热性的介质中的瞬态热分析受以下抛物线型偏微分方程的控制。

移动Kriging插值（MKI）[27][29]在无网格方法中受到青睐，因为它可以提供一个能够准确通过节点值并满足Kronecker δ性质的形状函数，并在整个问题域中表现出高阶连续性。在计算点X，温度场uhX通过MKI [27]由局部支持域中的n个节点的温度值近似表示：uhX = pTX_a + rTX_b，其中p(X) = [x, y, x_2, ...], T是多项式基向量；a = a_1, a_2, ..., a_m，T是多项式。数值例子：本节报告的所有时间测量都是在配备Core i7-12700KF处理器和32 GB RAM的标准台式计算机上获得的。使用的编程软件是MATLAB 2025b。为了减轻CPU调度和内存访问延迟等瞬时系统级波动的影响，每个节点的搜索时间独立测量五次并记录算术平均值。kd树搜索和全局搜索的计算仅在示例1中有所不同。

FGM的方形薄板：计算了一个FGM的方形薄板。方形板边长设为L=1，初始条件u_0、边界条件u_1、q_1和节点布置分别如图4和图5所示。材料密度ρ=1，导热系数k(X)和比热容c(X)的表达式如下：k_X = k_0 + βx^2，c_X = c_0 + βx^2，其中β是一个异质参数。k_0、c_0是原点处已知的导热系数和比热容值。

异质FGM中带孔的方形板：计算域是结构的四分之一（见图10）。方形板边长为L=4，中心圆形孔的半径为r=1。比热和密度的乘积为ρc=1。导热系数分别为k_11=x+1, k_22=y+1, k_21=k_12=0.5。图11显示了节点布局。初始条件u_0=0，边界条件、源项和解析解如下[46]。

示例6：三维不规则球体：在最后一个例子中，我们将MKI扩展到三维不规则体的瞬态热分析（见图30）。计算区域为[50]：?Ω = (x, y, z) / (x ≤ ηsinθsinφ, y ≤ ηsinθcosφ, z ≤ ηcosθ, 0 ≤ θ ≤ π, 0 ≤ φ ≤ 2π)，其中η=1+1/6sin6θ。导热系数为k=exp(z)，比热容为c=exp(z)，解析解的表达式如下：ux, y, z, t = sin(πx^2 + sin(πy^2 + sin(πz^2) + cos(πx^2 + cos(πy^2)) / (exp(-1/4βπ^2)t + 1)，x, y, z ∈ Ω, t ∈ [0, 3)。其余条件在文中推导。

结论：本文通过理论推导和六个基准数值例子，建立了一种结合移动Kriging插值（MKI）和kd树的高精度Galerkin无网格方法，用于FGMs的瞬态热分析问题。该方法构建了Galerkin弱形式时间步进方案，并实现了基于kd树的自适应邻域搜索策略。虽然保留了无网格方法的几何灵活性，但kd树用于节点搜索减少了计算成本。

作者贡献声明：吴希荣：撰写-审稿与编辑，撰写-初稿，可视化，验证，软件，资源，方法论，研究，形式分析，数据管理。凌静：撰写-审稿与编辑，撰写-初稿，验证，监督，资源，项目管理，方法论，研究，概念化。王晓斌：撰写-审稿与编辑，撰写-初稿，可视化，验证，监督，软件，方法论，形式分析。作者声明他们没有已知的竞争性财务利益或个人关系可能影响本文所报告的工作。

致谢：本研究得到了国家自然科学基金（项目编号11762005）的支持。