如何在有限元模拟中使用数值积分和高斯点

在有限元模型中，你可能会在多种情境下遇到数值积分和高斯点的概念。在本文中，我们将讨论在什么情况下，以及为什么使用数值积分。此外，还强调了在 COMSOL Multiphysics® 软件中检查和修改数值积分方案的可能性。最后，对高斯点自由度的使用进行了说明。

什么是数值积分？

在计算一般域上非平凡函数的积分时，我们必须采用数值方法。数值积分也称为数值求积，其本质是用求和代替积分，其中被积函数在多个离散点被采样，可以描述为

其中 xi 是积分点的位置，wi 是相应的权重因子。积分点通常称为高斯点，但是严格来说，这种命名法仅适用于高斯求积 法定义的积分点。在 COMSOL Multiphysics 中，真正的高斯求积用于一维，二维中的四边形单元，三维中的六面体单元等积分。其他情况下，需要使用其他类似的方案。

高斯求积

在高斯求积算法中，需要选择积分点的位置及其权重，以便精确地对阶次尽可能高的多项式进行积分。由于 N 次多项式包含 N + 1 个系数，而具有 M 个点的高斯点规则包含 2M 个参数（位置+权重），因此可以精确积分的多项式的最高阶次是   N = 2M-1。

在 COMSOL Multiphysics 中，积分阶次由可以精确积分的多项式的阶次指定，仅使用偶数。对于真正的高斯点积分，如上所示，精度始终是奇数次幂。软件指定为“4”的积分阶次实际上在这些单元形函数上使用精度为 5 阶的高斯求积。

高斯求积示例

要积分的函数。

4×4 积分点的情况。橙色表面是函数值为 1 的区域，绿色高斯点是对积分值有贡献的点。

从表中可以看出，计算这种不连续函数的精确积分是非常重要的。纯属偶然，2×2 积分方案给出了精确的答案，但收敛完全不单调。

为什么这很重要？在有限元分析中，你可能会遇到呈现明显局部梯度的场，例如相变问题或固体力学中塑性开始时的问题。在包含这种突变的单元上计算的积分可能具有显著的离散化误差。此外，解的收敛性也会减弱。当单个高斯点改变其状态时，解的微小变化可以显著改变计算残差。

在这种情况下，选择一个比用于将场离散化的形函数默认值更低的多项式阶次可能更好。较低的分辨率可以通过使用更密集的网格来补偿，其结果是不可避免的突变将被局限在具有较少积分点的较小单元。

此外，如果你在后处理过程中计算不连续函数的积分，要记住，数值积分收敛可能很缓慢。

弱贡献的积分

这里，波浪号（〜）表示虚变化。在 COMSOL Multiphysics 的表达式中，它由 test() 运算符表示。这里包含体积力 f 和边界牵引力 t，但也可能有其他类型的贡献。左侧将贡献于刚度矩阵，而右侧将贡献于载荷矢量（假设力与位移无关）。

检查稳态情况下 固体力学接口中 线弹性材料的弱表达式。

在上图中，你可以看到积分阶次的文本框，本例中，值为 4。COMSOL Multiphysics 中描述积分阶次的数字是可以精确积分的多项式的最高阶次，默认的积分阶次基于用于描述场（本例中为位移）的形函数的阶次。这里使用了默认的形函数阶次——二次，因此应力和应变在单元上基本上呈线性变化。所以，应力和应变变化的乘积是二次的，这表明 4 阶可能超过必要程度。为什么选择这个值将在下文讨论。

边界载荷的弱表达式。

这里，积分阶次也是 4。由于位移形函数是二次多项式，这意味着应该可以精确地积分牵引力不超过二次变化的载荷贡献，因此牵引力与位移变化的乘积是 4 阶。

关于弱表达式的一些微妙之处

雅可比通常是一个有理函数（分子和分母都是多项式），所以这种数值求积法甚至可能无法对它进行精确积分。因此，明智的做法是在选定的积分阶次中留出一些余量。顺便提一下，雅可比效应是严重失真的单元比具有理想形状的单元表现更差的原因之一。这篇关于 如何检查 COMSOL Multiphysics 中的网格 的博客文章包含了更多关于网格质量的信息。

即使单元形函数被称为“二次函数”，它也可能（在某些情况下）包含更高阶的项。用户界面中指示的形函数阶次显示了形函数中最高阶的完整多项式。使用比乍看起来更精确的积分规则的另一个原因是，多项式中可能存在一些高阶项。

对于某个贡献，你可能会在方程视图 中看到比预期更高的积分阶次的另一个原因是，可能需要确保刚度矩阵中的对称性。

修改积分阶次

你可以通过编辑文本框来更改方程视图 中任何弱表达式的积分阶次。你想要这样做有两个主要原因。

第一个也是比较明显的一个原因：你想提高精度。假设载荷（一般来说可以是力、热通量、电流等）与单元大小相比变化很快。那么，增加数值积分的阶次将会提高进入域的总力或通量的精度。但是，靠近应用边界条件的位置的局部解仍然不是很理想，原因是，它永远不会比单元的形函数所代表的更好。

然而，还有另一个有趣的情况：降阶积分，这意味着由于某种原因，积分阶次低于正式需要的阶次。其中一个原因是加快计算速度。在求解有限元问题时，大部分 CPU 时间花费在两个任务上：形成单元矩阵（装配）和求解大型线性方程组。求解方程所花费的时间比模型大小增加得更快，通常大约是单元数量的平方。在装配上花费的时间与模型大小成正比（实际上，与单元数和每个单元积分点数的乘积成正比）。对于非常大的模型，方程求解总是占主导地位，但对于在每个积分点有大量计算的中型非线性模型，可能值得考虑使用降阶积分。

轴对称 壳接口的虚功贡献。

积分耦合算子

添加积分算子。

添加积分算子时，需要进行三个主要选择：

1. 应采用积分的域、边界或边。

2. 积分阶次——这也为你提供了一个以精度换取速度的选项。请注意，实际的被积函数不仅是你提供的表达式，还要乘以从理想单元形状到真实单元形状转换的雅可比。

3. 框架——此项只有当存在不同框架（如动网格、变形几何和结构力学中的几何非线性）时才重要。简而言之就是：积分应该取变形的还是未变形的几何结构？

后处理过程中的积分

在结果评估期间添加积分节点。

高斯点形函数

为用户定义的因变量选择形函数类型。

当你想要存储局部状态时，高斯点数据 类型很有用。比如，在需要“记忆”的历史相关的非线性本构模型中就有这种情况。本构模型主要是在计算刚度矩阵和残差时访问的，因此在求解过程中，单元中唯一实际评估本构模型的位置是积分点。所以，将这种类型的数据准确存储在那里是有意义的。

内部使用高斯点数据的一个例子是在材料模型中存储非弹性应变，例如结构力学中的塑性和蠕变。

如果要添加高斯点变量以与内置物理场接口一起使用，通常应选择与用于计算相关弱表达式的积分阶次相同的积分阶次。如果不匹配，则必须在单元中的不同位置之间传输值，这种情况下精度和性能会降低。

gpeval 算子

如果变量存储在高斯点中，无论是内置还是由你定义，都有可能需要在单元上插入这些变量。这在后处理过程中尤为重要。默认情况下，当在单元中的另一个位置求值时，高斯点变量的值只是从最近的高斯点选取的。gpeval() 算子可用于将离散高斯点数据映射到连续场。在其最简单的形式中，算子以 gpeval(gporder, expression) 形式被引用，例如，gpeval(4,solid.epe)。有关更多详细信息，请参阅 COMSOL Multiphysics User’s Guide。

将 x 坐标作为高斯点数据存储。

弱贡献 (myX-nojac(X))*test(myX) 只表示：“将变量 myX 设置为当前值 X。”添加算子 nojac() 是为了防止 myX 和 X 之间出现双向耦合。当你只想为因变量赋值时，使用它是一个很好的做法。

高斯点变量（底部）和外推高斯点变量（顶部）图。箭头指示在结果评估期间高斯点数据默认是如何移动到拐角的。

本文内容来自 COMSOL 博客，点击“阅读原文”，阅读更多延伸文章。