有限元中的积分点结果是如何外推到节点的？

在有限元后处理中，我们经常会看到节点上的应力、应变、塑性应变等云图。例如 Abaqus 中可以直接显示节点上的 Mises 应力，看起来似乎每个节点天然就有一个应力值。但实际上，对于大多数基于 Gauss 积分的有限元，应力、应变以及大量材料状态变量并不是直接在节点上计算的，而是在积分点（integration point / Gauss point）上计算的。那为什么需要节点值呢？节点结果更适合工程后处理、结果展示和位置查询。后处理为了绘制光滑、连续的 contour，通常需要做以下操作：

\[\boxed{ \text{Integration Point} \rightarrow \text{Extrapolation} \rightarrow \text{Element Nodal} \rightarrow \text{Averaging} \rightarrow \text{Displayed Nodal Result} } \]

本文主要讨论其中最核心的一步：

如何把单元积分点上的变量外推（extrapolate）到单元节点上。

1. 问题定义

假设积分点上的某个标量变量为

\[\mathbf u_I = \begin{bmatrix} u_{I_1}\\ u_{I_2}\\ \vdots\\ u_{I_m} \end{bmatrix} \]

节点上的变量为

\[\mathbf u_N = \begin{bmatrix} u_{N_1}\\ u_{N_2}\\ \vdots\\ u_{N_n} \end{bmatrix} \]

我们希望找到一个线性变换

\[\boxed{ \mathbf u_N = \mathbf A\mathbf u_I } \]

其中 \(\mathbf A\in\mathbb R^{n\times m}\)。

就是所谓的积分点到节点的外推矩阵（extrapolation matrix）。问题于是变成：

这个矩阵 \(\mathbf A\) 是怎么得到的？

2. 从节点到积分点：形函数插值

先考虑一个更熟悉的问题。如果已知节点值 \(\mathbf u_N\)，如何计算积分点上的值？这正是有限元形函数最擅长做的事情。

对于任意单元内部位置 \(\boldsymbol\xi\)：\(u(\boldsymbol\xi) = \sum_{j=1}^{n} N_j(\boldsymbol\xi)u_{N_j}.\)

在第 \(i\) 个积分点 \(\boldsymbol\xi_i\) 上：\(u_{I_i} = \sum_{j=1}^{n} N_j(\boldsymbol\xi_i)u_{N_j}.\)

把所有积分点写在一起：

\[\begin{bmatrix} u_{I_1}\\ u_{I_2}\\ \vdots\\ u_{I_m} \end{bmatrix} = \begin{bmatrix} N_1(\boldsymbol\xi_1)&\cdots&N_n(\boldsymbol\xi_1)\\ N_1(\boldsymbol\xi_2)&\cdots&N_n(\boldsymbol\xi_2)\\ \vdots&&\vdots\\ N_1(\boldsymbol\xi_m)&\cdots&N_n(\boldsymbol\xi_m) \end{bmatrix} \begin{bmatrix} u_{N_1}\\ u_{N_2}\\ \vdots\\ u_{N_n} \end{bmatrix} \]

定义 \(N_{ij} = N_j(\boldsymbol\xi_i),\) 那么可以简写为 \(\boxed{ \mathbf u_I=\mathbf N\mathbf u_N }\)，其中 \(\mathbf N\in\mathbb R^{m\times n}\)。

3. 从积分点反推节点：外推

如果积分点数和节点数相同：\(m=n,\) 并且矩阵 \(\mathbf N\) 可逆，那么直接有

\[\boxed{ \mathbf u_N = \mathbf N^{-1}\mathbf u_I } \]

因此：

\[\boxed{ \mathbf A=\mathbf N^{-1} } \]

这就是最经典的积分点到节点外推矩阵。对于固定类型的单元，\(\mathbf A\) 只依赖参考单元，因此通常可以预先计算，而不需要对网格中的每一个单元重新求逆。数学上不能仅由 “积分点数 = 节点数” 推出 \(\mathbf N\) 一定可逆，但是商业有限元里满足 \(n_{\mathrm{GP}}=n_{\mathrm{node}}\) 的典型组合恰好都是低阶全积分单元，所以这种情况下 \(\mathbf N\) 都是可逆的。

4. 一个最简单的🌰：二节点一维单元

对于一维二节点线性单元：

\[N_1(\xi)=\frac{1-\xi}{2}, \qquad N_2(\xi)=\frac{1+\xi}{2} \]

两个 Gauss 点为 \(\xi_1=-\frac1{\sqrt3}, \qquad \xi_2=+\frac1{\sqrt3}\)。

因此：

\[\mathbf N= \begin{bmatrix} N_1(-1/\sqrt3)&N_2(-1/\sqrt3)\\ N_1(+1/\sqrt3)&N_2(+1/\sqrt3) \end{bmatrix} \]

代入形函数：

\[\mathbf N = \frac12 \begin{bmatrix} 1+\frac1{\sqrt3}&1-\frac1{\sqrt3}\\ 1-\frac1{\sqrt3}&1+\frac1{\sqrt3} \end{bmatrix} \]

求逆得到：

\[\mathbf A = \mathbf N^{-1} = \frac12 \begin{bmatrix} 1+\sqrt3&1-\sqrt3\\ 1-\sqrt3&1+\sqrt3 \end{bmatrix} \]

注意这里出现了一个非常有意思的现象：\(\frac{1-\sqrt3}{2}<0\)。也就是说，外推系数可以是负数，也可以大于 1。

例如节点 1：\(u_{N_1} = \frac{1+\sqrt3}{2}, \, u_{I_1} + \frac{1-\sqrt3}{2}u_{I_2}\)。这和普通“加权平均”完全不同。假设 \(u_{I_1}=100, u_{I_2}=200\)，那么 \(u_{N_1} \approx63.4\)，而 \(u_{N_2} \approx236.6\)。虽然积分点结果全部位于 \([100,200]\) 之间，但外推以后节点值却变成了 \(63.4,\,236.6\)。因此，节点外推值完全可能超出积分点结果的最大值和最小值。这对于理解商业有限元软件中的应力云图非常重要。

5. 当积分点数和节点数不相等怎么办？

前面的推导有一个很强的条件：\(n_{\mathrm{GP}}=n_{\mathrm{node}}\)。但实际有限元中，这个条件经常不成立。举个最简答的例子，C3D8R 八节点六面体减缩积分单元，它只有中心的一个积分点，却需要恢复 8 个节点上的值。这时 \(\mathbf N\) 明显不是方阵，因此不能简单写成：\(\mathbf A=\mathbf N^{-1}\)。这时候就必须回答一个更本质的问题：我们究竟假设积分点变量在单元内部服从什么样的场？

现在我们重新思考一下积分点外推的问题。它其实可以这么定义，假设单元内部存在一个场：

\[u_h(\boldsymbol\xi) = \sum_{j=1}^{n} N_j(\boldsymbol\xi)u_j \]

我们现在知道这个场在若干积分点上的值：

\[u_h(\boldsymbol\xi_1), u_h(\boldsymbol\xi_2), \ldots, u_h(\boldsymbol\xi_m) \]

所谓积分点外推，本质上就是：

利用这些离散采样值，先重构一个单元内部的场，然后在节点位置评价这个场。

所以可以把整个过程理解成：

\[\boxed{ \text{Gauss-point samples} \rightarrow \text{reconstructed field} \rightarrow \text{nodal evaluation} } \]

这样一来，为什么会有唯一性问题也就很好理解了。如果采样点提供的信息足够多，就能够唯一确定这个场；如果信息不够，就存在多个不同的场经过这些积分点时具有完全相同的值。这时候，所谓“外推到节点”就不再存在唯一答案，要唯一确定就必须额外引入某种假设。

6. Element Nodal \(\neq\) Unique Nodal

考虑两个相邻单元：

              shared node
                  ●
                 / \
                /   \
        element 1   element 2

虽然它们共享同一个几何节点，但是两个单元内部计算出来的应力场一般并不连续。因为对于普通 displacement-based (C^0)​ 有限元，我们主要保证的是位移场在单元之间连续，而并不要求应力场跨单元连续。所以：

\[\boxed{ \text{同一个几何节点可以对应多个 element nodal stress} } \]

这也就是为什么本文最开始的那个后处理流程里有个 averaging 的过程：我们需要把多个 element nodal values 合并成一个 unique nodal value，从而得到更加连续的 contour。

在 Abaqus 中，它还额外引入了 averaging threshold ，来控制共享节点处来自不同单元的 element-nodal values 是否进行平均。如果通过 averaging threshold 检查，Abaqus 才会把它们平均，如果没有通过 threshold，则不会把这两个 contribution 合并为同一个节点值，contour 在这里可以保持不连续。Abaqus 默认是 75%，averaging threshold 值越大，contour 越光滑。

7. 总结

最后，总结一下。

积分点到节点的外推，本质上是根据积分点处的离散采样值重构单元内部的场，并在节点位置对该场进行评价。当积分点数与节点数相同且 \(\mathbf N\) 可逆时，可以直接得到

\[\boxed{\mathbf u_N=\mathbf N^{-1}\mathbf u_I} \]

但当积分点数与节点数不相等时，问题不再是简单的矩阵求逆，而需要引入额外的场重构假设。因此，一般意义上的外推矩阵 \(\mathbf A\) 并不一定等于 \(\mathbf N^{-1}\)。

同时需要注意，Element Nodal 和 Unique Nodal 是两个不同层次的结果：前者由积分点结果在单元内部外推得到，后者还需要对相邻单元在公共节点上的结果进行 averaging。

除了本文介绍的基于形函数的外推之外，还有一种常见思路是最近节点映射（nearest-node mapping）：根据积分点与节点的空间对应关系直接进行结果映射。它与形函数外推的思路有明显区别，后续会单独写一篇介绍。