一种在无网格方法框架下用于梯度材料模型的计算高效数值积分方案

《Mathematics and Computers in Simulation》：A computationally efficient numerical integration scheme in the meshless method framework for graded material models

　　S.H. 斯里尼瓦桑|塔马赖·塞尔万·瓦苏|P.V. 杰亚卡蒂凯扬|赛·纳加·基肖尔·武特拉
印度瓦兰加尔国家技术学院机械工程系

**摘要**
本文将一种计算效率较高的数值积分方法——元素中点法（Element Midpoint Method）集成到径向点插值法（Radial Point Interpolation Method）框架中，用以研究在各种问题类型中处理功能梯度材料（Functionally Graded Materials, FGM）结构的计算效率。通过将结果与文献中报道的标准高斯求积法（Gauss Quadrature，阶数分别为3×3和6×6，GQ3和GQ6）进行比较来验证这些结果。此外，还记录了达到类似精度所需的计算时间，并将其与GQ3和GQ6的计算时间进行了对比。研究实例包括三种类型的问题：梁弯曲、自由振动和平面载荷下的板弯曲。为了验证结果，从文献中选取了线性、指数型和幂律型等材料梯度分布。研究结果表明，元素中点法在计算效率方面可以替代高斯求积法，同时具备较高的精度。

**引言**
材料研究已从均匀各向同性金属发展到合金、涂层、夹层结构、复合材料以及功能梯度材料（FGMs）等领域。在这些材料中，材料的性能根据几何形状（一维、二维和三维）进行定制，以实现稳健的设计。FGMs因其材料性能在整个结构中的连续变化而成为层状复合材料的先进替代品。最初，FGMs是为航空航天热障应用开发的，以承受高达1000°C的温度梯度[1][2]。与传统复合材料中的尖锐边界不同，FGMs能够实现材料的平滑过渡，从而最小化应力集中，提高抗热性和抗裂性，并减少层间脱离和横向剪切等界面问题[1][3][4]。尽管FGMs的连续材料梯度给精确建模和分析带来了挑战，但其独特的性能使其在航空航天、汽车和生物医学领域具有广泛应用潜力[1]。

另一方面，对于这类材料的建模和求解程序，需要考虑不同的几何形状、边界条件和载荷情况，同时保持足够的精度并减少计算时间。虽然解析方法可以提供连续且精确的解，但它仅适用于简单的几何形状，对于FGMs等复杂材料模型则显得力不从心。因此，引入了有限元法（Finite Element Method, FEM）、边界元法（Boundary Element Method, BEM）、等几何分析（Iso Geometric Analysis, IGA）和无网格方法（meshless methods）等数值方法。在这些数值方法中，高斯求积法是研究人员评估刚度矩阵的首选。尽管高斯求积法与FEM和BEM配合使用效果良好，但这主要得益于其所使用的多项式形状函数。然而，基于问题物理特性、解响应性质和相关奇异性的无网格方法需要高阶数值积分，从而导致更高的计算时间。

有限元法已被广泛用于分析功能梯度材料在各种载荷条件下的行为；不过此处仅简要回顾了几项代表性研究。Alshorbagy等人[5]研究了FGM梁的自由振动响应，Nguyen等人[6]采用全拉格朗日公式和牛顿-拉夫森迭代法对双相FGM夹层梁的大变形进行了非线性FEM分析，探讨了不同均质化方案的影响。Bian等人[7]开发了一个FEM框架，用于研究在不同边界条件和非局部参数下的FGM梁的静态弯曲、屈曲和自由振动。Boumezbeur等人[8]使用三维FEM模型和降阶积分方法研究了均匀载荷下FGM梁的静态响应，表明FEM能够有效捕捉到梁理论中忽略的三维效应。Benmalek等人[9]采用了混合等参FEM和降阶积分方法来准确高效地分析FGM梁的静态弯曲行为。Moita等人[10]基于圆锥台元素和傅里叶级数展开开发了FEM公式，用于轴对称FGM板壳结构的自由振动分析，通过一点高斯积分和惩罚技术准确捕捉薄板和厚板的横向剪切效应，同时降低了计算成本。Yasin等人[11]提出了一个两节点FEM，用于带有压电执行器的FG梁的动态分析和主动振动控制，结合了分层力学和Mori-Tanaka均质化方法，并采用六点高斯积分进行刚度和惯性评估。Mesbah等人[12]引入了准三维高阶FEM，用于厚FGM梁的自由振动和屈曲分析，通过选择性降阶积分抑制剪切锁定现象，同时保持了精度和收敛性。Burlayenko等人[13]在ABAQUS中使用了一维梁和三维固体模型进行了基于FEM的自由振动研究，表明全高斯积分的三维模型在更高阶模式和锥形几何形状下具有更高的精度，而一维模型在计算效率上仍具有优势。Martínez-Pa?eda[14]研究了基于高斯点和节点的材料梯度策略在FGM板平面载荷分析中的应用，比较了不同元素公式下的精度。研究显示，高阶高斯求积法（特别是全积分二次元素Q8）显著提高了精度，而基于节点的梯度方法引入了插值误差；降阶积分元素（Q4R、Q8R）虽然应力精度较低，但对材料梯度的敏感性较低。Vidal等人[15]提出了一种变量分离FEM框架，通过将问题分解为耦合的二维和一维线性子问题，使用3×3高斯积分进行元素矩阵评估，从而高效地进行了FGM板的非线性分析。

已经开发了多种无网格Galerkin公式用于在线性和非线性条件下分析功能梯度材料（FGMs）。Basant等人[16]结合Reddy的三阶剪切变形理论和Eringen的非局部模型，使用了无网格自然邻居Galerkin方法（NNGM）研究FG板的非线性弯曲响应，发现非局部参数引起的软化效应，并与FEM结果吻合良好。Liu等人[17]进一步利用再生核粒子方法（Reproducing Kernel Particle Method, RKPM）和全拉格朗日公式研究了几何非线性FGMs，通过高斯点评估空间变化的材料属性，从而提高了精度。Liu等人[18]引入了径向基再生核粒子方法（RRKPM），有效消除了与核相关的误差。Wang等人[19]基于移动最小二乘法和再生核形状函数开发了无网格Galerkin方法，用于六阶应变梯度梁和板问题，证明了一致的积分方案和平滑的节点导数在精度和效率上优于标准高斯积分。Thai等人[20]采用移动克里金法（Moving Kriging, MK）和非局部应变梯度理论研究FG纳米板的弯曲和自由振动，而Trai等人[21]扩展了MK框架，结合改进的FSDT方法对Kerr基础上的FG夹层板进行了振动分析。Vu等人[22]结合精细剪切变形理论和增强的MK插值方法，提高了弯曲、振动和屈曲分析的精度。

基于径向点插值的无网格方法因其Kronecker delta特性而被广泛应用于FGM问题，这种特性允许直接且精确地施加边界条件，而基于MLS或RK的方法则需要惩罚技术、拉格朗日乘数或变换方法。Van Do和Lee[23]提出了改进的RPIM用于非线性弯曲分析，Li和Liu[24]开发了无元素平滑RPIM以克服二维和三维固体力学中的背景网格依赖性。RPIM的应用包括功能梯度碳纳米管（FG-CNT）增强壳体的自由振动分析[25]、使用应变梯度理论的纳米板静态弯曲[26]，以及通过离散剪切间隙技术减轻剪切锁定的FG板非线性分析[27]。Vaghefi[28]提出了用于三维弹塑性后屈曲分析的Tchebychev-RPIM公式，强调了高斯求积阶数对收敛性和精度的影响。

在FEM、BEM、IGA和无网格公式等数值方法中，求解控制方程需要准确评估刚度和质量矩阵以及力向量，因此数值积分是一个关键环节。计算结果的精度很大程度上取决于所使用的积分方案，特别是在涉及材料属性连续变化的问题中，如FGMs[29]。高斯求积（GQ）是这些方法中最常用的积分技术，因其具有高精度、最佳收敛性和对平滑多项式形状函数的有效性，已在标准FEM公式中得到广泛应用[2][29]。尽管高斯求积具有鲁棒性，但在FGMs和无网格应用中存在局限性。精确的刚度和力评估通常需要高阶求积和大量积分点，这会增加浮点运算量和计算成本[2]。在无网格方法中，由于形状函数及其导数通常为非多项式和有理函数，高斯求积可能需要更多的积分点或无法满足积分一致性，从而导致精度降低[30]。此外，在材料不连续或存在裂纹的情况下，高斯求积因难以处理非光滑被积函数而效率低下且容易出错[2][4][29][31]。

为了克服高斯求积的这些局限性，研究人员探索了多种替代积分方法以提高精度和计算效率。Belytschko等人[32]提出了一种直接在节点处进行积分的节点积分方法。但后续研究表明，纯节点积分可能导致数值不稳定和刚度系统低估。为了解决这些问题[33]，他们引入了应变平滑策略并开发了稳定协调节点积分（Stabilized Conforming Nodal Integration, SCNI）方案，显著提高了节点积分的精度和稳定性。最近[34]将节点积分集成到RPIM框架（NI-RPIM）中，用于解决二维问题，展示了增强的鲁棒性、可靠性和数值稳定性。Jeyakarthikeyan等人[29]提出了两种计算效率较高的替代积分方法——元素中点法（Element Midpoint, EM）和元素边缘法（Element Edge, EE），以克服高斯求积在高斯求积分析中的高计算成本和局限性。他们的方法使用基于整数的采样替代了浮点高斯点，并结合Richardson外推法，实现了与高斯求积相当的精度，同时将计算成本降低了近三分之一，并有效控制了沙漏模式。作者进一步证明了这些方法在生成FGMs中四面体和六面体元素的元素矩阵方面的计算效率[31]，应用基于Richardson外推法的积分方案解决了厚FGM板的位移问题[2]，并扩展到使用基于应变能量的复合FGM板分析[4]，以及二维和三维结构的几何非线性分析[35]。Nguyen等人[36]使用EM和EE方法进行了二维和三维结构的几何非线性分析。文献表明，EM和EE方法在有限元法框架内被用于解决多种问题。随着FEM的成功，Vutla等人[37][38]将EM和EE方法引入RPIM无网格框架，用于解决断裂问题并证明了其计算效率。Kishore等人[30]开发了一种具有10个积分点的高阶元素中点法（EM-10），解决了断裂问题并预测了各种试样中的裂纹路径。他们通过与六阶高斯求积（GQ6，36个积分点）进行比较，证明了EM-10的速度比GQ6快3倍以上，且解决方案精度不受影响，即使对于包含多项式、三角函数和径向基函数的被积函数也是如此。

尽管文献表明元素中点法（EM）在FEM框架中得到了验证，但这些验证仅限于多项式函数。在本研究中，EM方法被应用于非多项式形状函数和非均匀材料模型。虽然EM方法在RPIM框架中用于断裂建模，但研究仅限于均匀材料。然而，研究FGMs在材料非线性条件下的计算效率（如静态弯曲、自由振动和平面载荷）非常重要且具有时效性，目前文献中尚未对此进行报道。目前的工作可以高效地用于对航空航天应用中使用的功能梯度梁和板进行结构分析，例如热屏障、飞机面板和涡轮机部件。无网格方法与电磁方法相结合，将有助于汽车和机械系统中梯度结构的振动和稳定性分析，这些分析需要多次模拟。在本文中，采用元素中点（EM）方法在静态和自由振动条件下分析功能梯度材料（FGMs），并进一步扩展到在无网格框架内解决二维平面载荷问题。文章的结构如下：第2节介绍了FGMs建模的问题表述和EM方法；第3节讨论了各种静态、自由振动和二维平面载荷问题的数值结果；第4节总结了研究得出的主要结论。

**问题表述**
在没有体力作用且施加牵引力和位移边界条件的情况下，应力平衡可以表示为方程（1）[39]：
$$
\rhu?i + \sigma_{ij,j} = 0, \quad \text{在 } \Omega \text{ 内}; \quad \sigma_n \cdot u = t_0; \quad \Gamma_t, \quad u = u_0; \quad \Gamma_u
$$
其中 $\rho$ 是密度，$\sigma_{ij}$ 是应力张量，$\Omega$ 是问题域，$t_0$ 和 $u_0$ 分别是规定的牵引力和位移边界条件，$n^\cdot$ 是单位外法向量。方程（1）在 $\Omega$ 内的Galerkin弱形式可以表示为：
$$
\int_{\Omega} (L \delta u) T_d(Lu) d\Omega - \int_{\Gamma} \delta u T_t d\nu = 0
$$

**数值结果与讨论**
在本节中，提供了以下数值示例来验证基于第2节中给出的RPIM建模开发的Matlab代码，并展示EM方法的计算效率：
1. 图3所示的悬臂FGM梁在梁全长上受到均匀载荷（UDL）作用的静态分析；
2. 图5所示的简支FGM梁在梁全长上受到均匀载荷（UDL）作用的静态分析；

**结论**
二维功能梯度材料使用径向点插值方法进行建模，以解决固体力学问题，如梁弯曲、自由振动和平面载荷下的板。元素中点方法用于积分非多项式函数，作为高斯求积法的替代方案。通过对上述各种问题的研究，可以得出以下结论：
1. 元素中点方法无论对于哪种类型的网格都显示出良好的准确性。

**作者贡献声明**
S.H. Srinivasan：撰写——原始草稿，可视化，验证，软件，方法论，形式分析。
Thamarai Selvan Vasu：撰写——审阅与编辑，撰写——原始草稿，可视化，验证，监督，方法论，形式分析，概念化。
P.V. Jeyakarthikeyan：撰写——审阅与编辑，撰写——原始草稿，可视化，监督，研究，概念化。
Sai Naga Kishore Vutla：撰写——原始草稿，可视化，验证，软件，资金数据。

我们希望感谢国防研究与发展组织（DRDO）通过CARS项目为我在Warangal国家技术学院的研究提供的财政支持。

**致谢**
作者们衷心感谢印度理工学院马德拉斯分校航空航天工程系的前教授G. Subramanian教授提供的宝贵指导。