用于纤维增强复合材料中观尺度分析的高阶嵌入有限元方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：High-order embedded finite elements for the mesoscale analysis of fiber-reinforced composites

　　摘要：复合材料的多尺度特性要求在各个材料尺度上进行精确建模，以确保在不同层级间实现有效的信息传递。基于经典理论（如Reissner–Mindlin理论）的有限元模型通常被用于描述层合结构的宏观行为，而三维方法则多用于分析纤维-基体或微观层面的行为。然而，在介观尺度建模方面仍存在不足：经典模型精度不够，三维模型则需耗费大量计算资源。因此，高效的介观尺度分析工具至今仍未得到充分发展。为弥补这一缺陷，本研究提出了一种具有可调精度的介观尺度分析数值策略。该方法采用通过Carrera统一公式生成的、具有可变运动学的嵌入有限元，通过拉格朗日乘数将宿主域与嵌入域的运动学场相互耦合。此外，该方法还解决了现有商业嵌入有限元框架的两个局限性：体积冗余以及仅能耦合平移自由度的问题。实验结果与现有研究高度一致，证明了所提出的基于CUF的嵌入有限元方法的准确性与高效性。

引言：近年来，计算能力的提升，尤其是高性能计算集群的发展，大大降低了在多个结构尺度上对材料行为进行建模的难度。不过，对于大规模结构而言，直接数值模拟的计算成本仍然过高。这一问题在复合材料中尤为突出，因为这类材料需要根据具体应用需求，采用多种不同的增强相与基体组合。在宏观尺度上，人们通常采用经典层合板理论或一阶剪切变形理论来进行分析。而在介观和微观尺度上，若要实现精确建模，则需要采用三维模型来描述各组分间的应力梯度，尤其是在进行详细或渐进性失效分析时。目前，针对复合材料微观和介观行为建模的方法主要有解析法、半解析法以及数值法三类。这些方法适用于相对简单的代表性体积单元，如方形或六边形排列、夹杂物或随机分布的纤维结构。而对于纺织复合材料这类更为复杂的结构，对其纤维-基体代表性体积单元进行离散化处理会显著增加难度。基于体素的有限元方法因能大幅减少单元生成前的预处理时间而被广泛采用，但这类模型在界面处的应力场预测往往不够准确。为确保位移连续性，人们通常会构建所有组件在界面处共享节点的符合网格，但当需要高精度应力预测时，这种方式往往会引入过多的自由度。嵌入有限元方法作为一种有效的替代方案，将代表增强相的网格置于基体网格内部，并通过约束嵌入网格的自由度，使其值与宿主网格的插值结果保持一致。尽管嵌入有限元方法在过去十年中受到了更多关注，但其基本原理其实早已确立。MSP技术最初由Ortiz等人提出，旨在提升等参元素在模拟局部变形时的性能。应变局部化现象指的是由于材料的本构特性，变形集中在狭窄区域内，如剪切带中。为模拟这类局部区域，Ortiz等人引入了不在元素边界保持C0连续性的形状函数。Belytschko等人进一步发展了这一理念，将局部区域的应变场视为连续的，而非不连续的。Fish和Belytschko则将这种方法应用于大变形问题，推导出了局部元素中节点力与刚度矩阵的一致表达式。他们的研究表明，剪切带的取向会显著影响元素的刚度，而较小的剪切带，如颈缩现象中所观察到的那种，对应的刚度值更低。Fish随后提出了叠加有限元方法，即s版本，用于处理具有陡峭应力梯度的问题，如固体中的裂纹或流体中的冲击波。在这些区域，若采用均匀网格划分，则需要大量的元素，从而大幅增加计算成本。与在保持多项式阶数不变的同时细化网格的h方法，或在固定网格上提高多项式阶数的p方法不同，s版本会在应力梯度较高的区域，如局部失效区，叠加更高阶数的元素。后来，Fish和Markolefas进一步扩展了这一技术，用于计算多层层合结构中的应力状态，能够处理奇异解与非奇异解。同样，Fish和Guttal也将这一框架应用于层合壳结构，通过在全局网格上叠加局部的一组三维元素，并结合二维等效的单层理论。上述关于叠加方法的研究主要聚焦于复合层合结构的层级分析。然而，要全面理解复合材料中各相的行为，多尺度分析是不可或缺的。为此，人们开发出了能够直接模拟纤维与基体的高效介观尺度建模技术。例如，Jiang等人提出了域叠加技术，用于分析编织复合材料。在这种方法中，纤维与基体分别进行网格划分，再通过基于插值的方程将二者耦合起来形成整体模型。纤维与基体之间的位移通过形状函数联系起来，这些函数可将局部纤维坐标映射到全局基体网格上。实验结果显示，该模型在应力与应变分布方面与全有限元模型的结果高度一致，不过该研究仅限于较为简单的几何形状和拉伸载荷条件。Iarve等人提出了独立网格方法，该方法结合了直接网格划分与基于体素的网格划分技术，可用于模拟复杂的复合结构几何形态。在这种方法中，纤维与基体分别独立网格化，然后再进行叠加，通过将重叠区域的相应形状函数设为零来消除多余的体积部分。在过去的十年中，还有许多研究采用了嵌入方法。特别是Tabatabaei等人的研究，评估了嵌入方法在纤维增强复合材料介观尺度分析中的适用性及其与传统有限元方法的精度对比。他们还明确指出了基体相中的体积冗余问题，并在计算刚度矩阵及结果后处理过程中提出了相应的解决方法。Liu等人和Tabatabaei等人后来进一步发展了这一方法，分别用于考虑纤维增强复合材料和普通编织复合材料代表性体积单元中的界面损伤与基体损伤等非线性现象。这些研究表明，域叠加技术不仅有助于简化模型构建并提升网格质量，还能使分析结果与传统有限元分析及实验测试的结果保持一致。尽管嵌入方法在网格划分方面具有灵活性，但上述研究大多仍采用三维实体模型来精确描述应力分布。实体元素建模存在一些限制，比如对元素长宽比和偏斜度的要求较高，往往需要大量元素才能获得满意的分析结果。不过，近期有研究表明，一维有限元也能足够好地模拟复合材料中的纤维相，且其结果与实验值高度吻合。例如，Salviato等人证明了Timoshenko梁单元的适用性。这种方法基于Cusatis等人提出的格点离散粒子模型，尤其是该模型在纤维增强混凝土中的应用，其中用于离散化纤维的节点同时也作为表示基体相的面的节点。虽然[19]中的方法并非嵌入技术，但与之有一定相似之处。不过，其与前述方法以及本文提出的方法不同，其变分形式并未考虑体积冗余项。商业有限元软件中也提供了嵌入功能，但这些功能存在若干局限性：其一，它们仅能约束平移自由度，当将一维或二维元素嵌入三维实体宿主中时，旋转自由度则不受约束；其二，它们仅适用于位移场分析；其三，它们无法处理体积冗余问题。第一个局限性在岩土工程领域已有相关研究。Turello等人开发了一种嵌入梁单元，用于在理想粘附假设下模拟土壤与桩基之间的相互作用，通过拉格朗日乘数来实现运动学上的兼容性。后来，Granitzer等人改进了这一方法，去掉了Turello等人使用的投影算子，同时引入了惩罚势函数，从而更有效地实现梁单元与实体单元之间的耦合。这两项研究都采用了Timoshenko梁理论，且仅考虑了圆形截面桩基，因此适用范围有限。为解决这一问题，本文提出了一种新的框架，能够处理任意截面的梁结构，以及从低阶到高阶的各种运动学模型。此外，本文还展示了如何将表示纤维的一维有限元嵌入到作为粘结基体或实际结构的另一维有限元中，这两种结构具有不同的力学描述方式。为实现这一目标，本文采用了Carrera统一公式作为结构理论的生成工具。无论采用何种运动学模型，CUF都能以简洁的形式表达控制方程。此前已有研究利用CUF实现了不同运动学模型之间的耦合，例如Carrera等人通过点态拉格朗日乘数，将多种一维运动学模型耦合起来，用于分析紧凑型、薄壁结构以及类似桥梁结构的局部区域。Pagani等人和Zhang等人则通过拉格朗日乘数，将基于键合的三维peridynamics方法与高阶一维及二维有限元相结合，用于进行准静态断裂分析。本研究进一步发展了作者在[29]中提出的初始CUF-EFE公式，该公式在之前的研究中仅考虑了耦合过程中的纯平移自由度。本文提出了一种嵌入有限元框架，其中宿主域与嵌入域的运动学描述均通过CUF生成，并通过拉格朗日乘数实现耦合，无需对嵌入体的形状或所选结构理论做任何特殊假设。本文的结构如下：第二节介绍CUF框架及其在线性弹性问题中的应用；第三节详细阐述宿主域与嵌入域之间的运动学耦合机制；第四节展示不同纤维-基体组合的仿真结果；第五节总结主要研究结论。

片段摘录：

结构理论的统一表述：考虑图1所示的笛卡尔参考系，其中三维位移向量可表示为：u(x,y,z)={ux,uy,uz}T，而应变与应力张量则采用Voigt记号表示为：????????={?xx,?yy,?zz,2?xz,2?yz,2?xy}T，σ={σxx,σyy,σzz,σxz,σyz,σxy}T。在位移较小的情况下，位移与应变可通过微分算子D关联起来，该算子的表达式为：D=??x000??y000??z??z0??x0??z??y??x??y0。随后，应力与应变之间也可通过某种方式建立关联。

宿主域与嵌入域之间的运动学耦合：不同子域之间的连接方式早已被研究过，人们采用了拉格朗日乘数、惩罚函数以及Nitsche方法等技巧来解决弹性问题。不过，这些方法通常用于将使用相同运动学模型的不相交子域进行组合。而在本文中，我们则将多个物体视为嵌入在某个宿主物体中，因此需要对方程(10)中的PVD公式进行重新推导，从而允许使用不同的运动学模型。

单纤维单元：该方法的正确性首先由Tabatabaei等人[15]通过对受均匀横向位移作用的方阵排列单元进行验证而得到确认。此时纤维的直径为d=1，单元宽度为b=1.429，长度为L=10。纤维与基体成分的材料属性列于表1中。所采用的边界条件如下：在y=0和y=L的面处uy=0，在z=±b/2处uz=0，在x=?b/2处ux=0，在x=b/2处ux=b/10。图4展示了相关的图形结果。

结论：本研究提出了一种嵌入有限元框架，该框架能够结合多种运动学描述方式，用于在不同边界条件与载荷作用下模拟纤维增强结构。用于描述各组成成分的结构理论是通过Carrera统一公式直接推导得到的。增强相采用不同阶数的泰勒展开，而基体或均匀化层合层则同时采用泰勒展开与拉格朗日展开两种方法。作者贡献声明：Alberto Racionero Sánchez-Majano：形式分析、资金获取、研究工作、软件使用、初稿撰写、审稿与编辑。Alfonso Pagani：资金获取、资源提供、审稿与编辑。

利益冲突声明：作者声明不存在任何可能影响本文研究结果的已知财务利益或个人关系。

致谢：第一作者得到了UC3M自主研究计划中“青年博士研究活动资助”项目的支持。同时，也感谢意大利大学与研究部通过FARE-III项目（LOUD项目，编号R20EENHZEJ）提供的支持。

Alberto Racionero Sánchez-Majano | Alfonso Pagani