在有限元模型中，你可能会在多种情境下遇到数值积分和高斯点的概念。在本篇博客文章中，我们将讨论在什么情况下，以及为什么使用数值积分。此外，还强调了在 COMSOL Multiphysics® 软件中检查和修改数值积分方案的方法。最后，对高斯点自由度的使用进行了说明。
目录
什么是数值积分？
在计算一般域上非平凡函数的积分时，我们必须采用数值方法。数值积分也称为数值求积，其本质是用求和代替积分，其中被积函数在多个离散点被采样，可以描述为
其中xi是积分点的位置，wi是相应的权重因子。积分点通常称为高斯点，但是严格来说，这种命名法仅适用于高斯求积法定义的积分点。在 COMSOL Multiphysics 中，真正的高斯求积用于一维，二维中的四边形单元，三维中的六面体单元等积分。其他情况下，需要使用其他类似的方案。
高斯求积
在高斯求积算法中，需要选择积分点的位置及其权重，以便精确地对阶次尽可能高的多项式进行积分。由于 N 次多项式包含 N + 1 个系数，而具有 M 个点的高斯点规则包含 2M 个参数（位置+权重），因此可以精确积分的多项式的最高阶次是 N = 2M-1。
高斯求积对于可由一定程度的多项式很好地进行近似的积分场非常有效。一维中一阶高斯求积的积分点和权重如下表所示。积分取自归一化区间 [-1,1]。
| 阶次（M） | 精度（N） | 位置（xi） | 权重（wi） |
|---|---|---|---|
| 1 | 1 | 0 | 2 |
| 2 | 3 | ±0.577 | 1 |
| 3 | 5 | 0, ±0.775 | 0.889 (= 8/9), 0.556 (= 5/9) |
| 4 | 7 | ±0.340, ±0.861 | 0.652, 0.348 |
| 5 | 9 | 0, ±0.538, ±0.906 | 0.569, 0.479, 0.237 |
在 COMSOL Multiphysics 中，积分阶次由可以精确积分的多项式的阶次指定，仅使用偶数。对于真正的高斯点积分，如上所示，精度始终是奇数次幂。软件指定为“4”的积分阶次实际上在这些单元形函数上使用精度为 5 阶的高斯求积。
高斯求积示例
作为高斯求积的一个示例，设想一个函数
{2}
\right )
在 -1≤x≤1，-1≤y≤1 的正方形上，积分为 1。如下图所示，该函数在域上的分布相当复杂。
要积分的函数。
下表显示了不同高斯求积阶次的结果。只要函数可以用一定阶次的多项式合理地近似，积分的值就会快速收敛。
| 积分点 | 精度 | 积分阶次 | 值 | 注释 |
|---|---|---|---|---|
| 1 | 1 | 0（或 1） | 2.9958 | 仅使用质心处的值明显高估了积分 |
| 2×2 | 3 | 2（或 3） | 0 | 在 2×2 规则中，高斯点位于\pm \frac{1}{\sqrt{3}}，余弦函数为 0（运气太差了！） |
| 3×3 | 5 | 4（或 5） | 1.1519 | |
| 4×4 | 7 | 6（或 7） | 0.9887 | |
| 5×5 | 9 | 8（或 9） | 1.0005 | |
| 6×6 | 11 | 10（或 11） | 1.0000 | |
| 7×7 | 13 | 12（或 13） | 1.0000 | |
| 8×8 | 15 | 14（或 15） | 1.0000 |
然而，多项式的最优行为高斯求积有一个缺点：不能很好地对非常不连续函数进行积分。假设我们要对如下函数进行积分：在上述相同正方形上，y < -2x – 1 时，函数值为 f(x,y) = 1，其他情况下，f(x,y) = 0。由于函数在正方形四分之一区域上的值为 1，并且正方形面积为 4，我们可以立即看出精确积分应该为 1。结果如下表所示。
| 积分点 | 值 | 注释 |
|---|---|---|
| 1 | 0 | 只有（0,0）一个点上函数的计算结果为 0 |
| 2×2 | 1 | 四分之一点的计算结果为 1，权重为 1 |
| 3×3 | 0.8025 | 两个点的计算结果为 1，权重分别为\frac{25}{81}，\frac{40}{81} |
| 4×4 | 1.2269 | 五个点的计算结果为 1 |
| 5×5 | 1.0325 | |
| 6×6 | 1.0918 | |
| 7×7 | 0.9892 | |
| 8×8 | 0.9961 |
4×4 积分点的情况。橙色表面是函数值为 1 的区域，绿色高斯点是对积分值有贡献的点。
从表中可以看出，计算这种不连续函数的精确积分是非常重要的。纯属偶然，2×2 积分方案给出了精确的答案，但收敛完全不单调。
为什么这很重要？在有限元分析中，你可能会遇到呈现明显局部梯度的场，例如相变问题或固体力学中塑性开始时的问题。在包含这种突变的单元上计算的积分可能具有显著的离散化误差。此外，解的收敛性也会减弱。当单个高斯点改变其状态时，解的微小变化可以显著改变计算残差。
在这种情况下，选择一个比用于将场离散化的形函数默认值更低的多项式阶次可能更好。较低的分辨率可以通过使用更密集的网格来补偿，其结果是不可避免的突变将被局限在具有较少积分点的较小单元。
此外，如果你在后处理过程中计算不连续函数的积分，要记住，数值积分收敛可能很缓慢。
弱贡献的积分
在每个有限元内，有各种表达式需要进行积分，以便形成刚度矩阵、质量矩阵、载荷矢量和残差等。以固体力学为例，标准弱（或变分）公式对应于虚功原理：内应力作用于虚应变变化所做的虚功等于外力作用于相应的虚位移变化所做的虚功。
这里，波浪号（〜）表示虚变化。在 COMSOL Multiphysics 的表达式中，它由test()运算符表示，其中包含体积力f和边界牵引力t，也可能还包含其他类型的贡献。左侧将贡献于刚度矩阵，而右侧将贡献于载荷矢量（假设力与位移无关）。
为了检查 COMSOL Multiphysics 中有限元实现的公式，你需要启用方程视图。
打开方程视图。
如果我们在固体力学接口中查看材料模型（比如线弹性材料）下的方程视图，会看到一个名为弱表达式的栏，你可以从中看到用于形成各种矩阵的表达式，例如，刚度矩阵。
检查稳态情况下固体力学接口中线弹性材料的弱表达式。
在上图中，你可以看到积分阶次的文本框，本例中，值为 4。COMSOL Multiphysics 中描述积分阶次的数字是可以精确积分的多项式的最高阶次，默认的积分阶次基于用于描述场（本例中为位移）的形函数的阶次。这里使用了默认的形函数阶次——二次，因此应力和应变在单元上基本上呈线性变化。所以，应力和应变变化的乘积是二次的，这表明 4 阶可能超过必要程度。为什么选择这个值将在下文讨论。
第二个例子，我们来看固体力学中的边界载荷。
边界载荷的弱表达式。
这里，积分阶次也是 4。由于位移形函数是二次多项式，这意味着应该可以精确地积分牵引力不超过二次变化的载荷贡献，因此牵引力与位移变化的乘积是 4 阶。
关于弱表达式的一些微妙之处
只有当单元具有理想形状（例如，没有弯曲边界）时，上述讨论内容才完全正确。如果你深入研究这个理论，积分实际上也会包含局部比例因子（雅可比），它来自实际单元几何结构与名义单元几何结构之间的转换。例如，如果在二维四边形单元上进行积分，则在理想正方形 -1≤ξ≤1，-1≤η≤1 上进行数值计算，
\displaystyle \int_{\Omega_e} f(x,y)\;dA = \int_{-1}^{1} \int_{-1}^{1} f(\xi, \eta) \left | \dfrac{\partial \mathbf x}
{\partial \boldsymbol \xi}
\right | \; d\xi d\eta
雅可比通常是一个有理函数（分子和分母都是多项式），所以这种数值求积法甚至可能无法对它进行精确积分。因此，明智的做法是在选定的积分阶次中留出一些余量。顺便提一下，雅可比效应是严重失真的单元比具有理想形状的单元表现更差的原因之一。这篇关于在 COMSOL Multiphysics 中检查网格的博客文章包含了更多关于网格质量的信息。
即使单元形函数被称为“二次函数”，它也可能（在某些情况下）包含更高阶的项。用户界面中指示的形函数阶次显示了形函数中最高阶的完整多项式。因此，使用比看起来更精确的积分规则的另一个原因是，多项式中可能存在一些高阶项。
对于某个贡献，你可能会在方程视图中看到比预期更高的积分阶次的另一个原因是，可能需要确保刚度矩阵中的对称性。
修改积分阶次
你可以通过编辑文本框来更改方程视图中任何弱表达式的积分阶次，这样做有两个主要原因。
第一个也是比较明显的一个原因：提高精度。假设载荷（一般来说可以是力、热通量、电流等）与单元大小相比变化很快，那么增加数值积分的阶次将会提高进入域的总力或通量的精度。但是，靠近应用边界条件的位置的局部解仍然可能不是很理想，原因是，它永远不会比单元的形函数所代表的更好。
然而，还有另一个有趣的情况：降阶积分。这意味着由于某种原因，积分阶次低于正式需要的阶次，其中一个原因是加快计算速度。在求解有限元问题时，大部分 CPU 时间花费在两个任务上：形成单元矩阵（装配）和求解大型线性方程组。求解方程所花费的时间比模型大小增加得更快，通常大约是单元数量的平方。在装配上花费的时间与模型大小成正比（实际上，与单元数和每个单元积分点数的乘积成正比）。对于非常大的模型，方程求解总是占主导地位，但对于在每个积分点有大量计算的中型非线性模型，可能值得考虑使用降阶积分。
降阶积分也是一种数值装置，有时用于消除人工刚度，这种人工刚度可能出现在某些单元公式中，这种现象通常称为锁定。在二维轴对称的壳接口中可以找到为此目的使用降阶积分的例子。在虚功方程中，一些项采用阶次 2 进行积分，另一些项采用阶次 4 进行积分。如果在任何位置都使用完全积分，当壳厚度变小时，单元实际上会变得太硬。通过使用降阶积分，可以有意地丢弃应变能中一些有问题的项。
轴对称壳接口的虚功贡献。
编者注：以下部分是于2023年5月4日新增的内容，包括关于 COMSOL Multiphysics 6.0 版本新功能的信息。
内置的减缩积分
一些结构力学的接口包括可以直接从材料模型的设置中选择减缩积分选项。
线弹性材料设置窗口的正交设置部分。
与在方程视图中编辑积分阶次相比，这种方法有几个优点：
- 改变可以在一个地方进行，并自动传递到任何可能的子节点。
- 你不必担心什么是适当的减缩积分方案。
- 对于某些单元形函数，如果降低积分阶次，刚度矩阵将变得奇异，可以通过内置的沙漏稳定性功能进行调整。
- 有些材料模型（如塑性），在积分点存储了局部状态，这是自动与积分规则同步的。
积分耦合算子
在“模型开发器”的定义下，你可以创建积分算子，用来定义作为问题公式一部分的全局变量，但它们也可以在结果评估期间显式地用在表达式中。
添加积分算子。
添加积分算子时，需要进行三个主要选择：
- 应采用积分的域、边界或边。
- 积分阶次——这也为你提供了一个以精度换取速度的选项。请注意，实际的被积函数不仅是你提供的表达式，还要乘以从理想单元形状到真实单元形状转换的雅可比。
- 框架——此项只有当存在不同框架（如动网格、变形几何和结构力学中的几何非线性）时才重要。简而言之就是：积分应该取变形的还是未变形的几何结构？
积分算子的设置。
后处理过程中的积分
如果要在后处理期间计算积分，可以使用两个选项：使用积分算子（如上所述）或在派生值下添加积分节点。在很大程度上，选择随意。但是，在积分节点中，你无法明确选择积分框架；它是根据数据集节点中的框架选择推断出来的。
在结果评估期间添加积分节点。
数据集的设置，从中可以选择结果解释的框架。
如果框架选择很重要，你可能应该使用积分算子将出现细微错误的风险降至最低。如果在求解完成后添加积分算子，则必须先运行更新解，然后才能使用新算子。
更新解。
请记住选择足够高的积分阶次。特别是，如果要积分的表达式是强非线性或不连续的，更加要注意。不连续表达式的一个常见特例是布尔表达式。例如，如果在传热分析后，你想计算温度高于某个值的体积，你可以使用像T>1066[K]这样的被积函数，该表达式在满足条件的情况下计算结果为 1，在其他情况下计算结果为 0，因此对它进行积分将给出满足条件时的体积。但是，这两个值之间的边界通常会穿过单元。
积分精度提高后的布尔表达式积分。
高斯点形函数
如果要对多物理场模型进行扩展，有时需要添加用户定义的自由度（因变量）。执行这个操作时，你必须选择用于表示它们的形函数的类型，其中一个选项是高斯点数据。所有其他选项给出不同类型的场，这些场在单元上为连续分布，在相邻单元之间可以是连续的，也可以是不连续的。“形函数”的高斯点数据类型完全不同，它仅在每个高斯点存储一个值，但与单元中其他位置的值没有关联。
为用户定义的因变量选择形函数类型。
当你想要存储局部状态时，高斯点数据类型很有用。比如，在需要“记忆”的历史相关的非线性本构模型中就有这种情况。本构模型主要是在计算刚度矩阵和残差时访问的，因此在求解过程中，单元中唯一实际评估本构模型的位置是积分点。所以，将这种类型的数据准确存储在那里是有意义的。
内部使用高斯点数据的一个例子是在材料模型中存储非弹性应变，例如结构力学中的塑性和蠕变。
决定存储高斯点数据后，就需要选择单元阶次。在高斯点数据情况下，这与上面讨论的积分阶次相同。存储高斯点数据的成本（在内存和 CPU 时间方面）与一维中的选定阶次、二维中的平方和三维中的三次方成正比。
选择积分点模式。
如果要添加高斯点变量以与内置物理场接口一起使用，通常应选择与用于计算相关弱表达式的积分阶次相同的积分阶次。如果不匹配，则必须在单元中的不同位置之间传输值，这种情况下精度和性能会降低。
gpeval 算子
如果变量存储在高斯点中，无论是内置还是由你定义，都有可能需要在单元上插入这些变量。这在后处理过程中尤为重要。默认情况下，当在单元中的另一个位置求值时，高斯点变量的值只是从最近的高斯点选取的。gpeval()算子可用于将离散高斯点数据映射到连续场。在其最简单的形式中，算子以gpeval(gporder, expression)形式被引用，例如，gpeval(4,solid.epe)。有关更多详细信息，请参阅COMSOL Multiphysics User’s Guide。
为了进行说明，请设想以下示例：x坐标（范围从 0 到 3）作为高斯点数据存储在小型三单元模型中。这是通过添加辅助因变量来实现的，如下所示。
将x坐标作为高斯点数据存储。
弱贡献(myX-nojac(X))*test(myX)只表示：“将变量myX设置为当前值X。”添加算子nojac()是为了防止myX和X之间出现双向耦合。当你只想为因变量赋值时，使用它是一个很好的做法。
如果高斯点变量myX随后被绘制为曲面图（没有单元之间的平均值），结果将是不连续的。在每个单元中，高斯点变量被移动到最近的角，然后插值到单元上。相反，如果绘制表达式gpeval(2,myX)，我们将检索精确的x坐标分布。
高斯点变量（底部）和外推高斯点变量（顶部）图。箭头指示在结果评估期间高斯点数据默认是如何移动到拐角的。
后续操作
单击下面的按钮，了解有关 COMSOL® 软件中可用功能的更多信息：
评论 (1)
雄 纪
2021-11-24最后两点“高斯点形函数”和“gpeval 算子”看不懂啊