有限元分片试验 | 从理论到手搓代码step-by-step数值实现！

在我们自己有限元编程的过程中，数值精度验证方式通常与商软、试验、理论解做对比，在本次中给大家分享一个除了上述方法外的精度验证方法—分片试验（patch test）有时也称拼片试验。

第一次接触"分片试验"概念是在学习曾攀老师的《有限元分析基础教程》，当时只是有个印象，具体怎么展开还是不知如何下手，直到最近读到了张雄老师的《有限元法基础》，才梳理一下分片试验的数值实现过程，现在分享给大家，希望能帮助到有需要的小伙伴。

本次推文主要有以下方面内容：

1. 理论介绍
2. 如何分片试验
3. 数值算例
4. 数值实现
5. 云图解释

理论介绍

为保证有限元解的收敛性，离散的单元必须满足一致性条件和稳定条件。

• 一致性条件：当单元尺寸趋于0时，刚度方程$Kd=f$可以代表原微分方程和边界条件；
• 稳定性条件：要求刚度方程$Kd=f$的解是唯一的，即整体刚度矩阵$K$是非奇异的。

分片试验实质上是在测试单元的稳定性条件，对于以上概念若觉得抽象，可暂且不看，直接看如何数值实现！

试验通过的标准：若干个单元组成的任意分片都能够表示试验模型的常应变状态。

什么叫做常应变状态？

可参考常应变三角形单元章节内容，即求解域内的节点位移模式均为线性多项式：
$$u(x,y)=a+bx+cy$$

分片试验一般选取4-8个任意形状的单元组成分片，见下图，在求解域内的各个节点上按照上式施加位移或力，得到的结果和线性多项式给出的解需要在计算机精度内保持一致。

精度要求

双精度浮点数最多能表示16位有效数字，因此采用双精度浮点运算的分片试验计算结果和线性多项式给出的解之间的误差必须在$10^{-15}-10^{-16}$之间，若不满足此精度要求，则需要检验程序的正确性。

值得注意的是：分片试验的单元的几何形状必须为非规则的，因为有时单元可以通过规则的分片试验，但不能通过非规则的分片试验！

接下来重点来了！编程中理解以下概念即可。

如何分片试验

分片试验有3种可能实施的方式：

• 试验A：在求解域内的所有节点施加按照线性多项式给定位移，检验内部节点是否满足平衡条件，是有限元解收敛的必要条件；
• 试验B：在求解域内的边界节点按照线性多项式施加给定位移，检验由有限元计算的节点位移是否与按照线性多项式确定的位移相同，是有限元解收敛的必要条件；
• 试验C：在求解域内的边界节点只施加能消除分片刚体位移的最少量本质边界条件（左侧边界节点固支，右侧约束y方向的自由度），其余边界节点按照线性多项式确定的面力施加自然边界条件，检验由有限元计算的节点位移是否与按照线性多项式确定的位移相同，是有限元解收敛的充分必要条件；

在进行分片试验时，需要给定线性多项式中的多项式系数，这些系数的取法有无穷中可能，可根据需要进行选项，如最简单的单向拉伸应力状态。

数值算例

现使用5个平面应力四节点等参单元组成一个分片试验，如下图所示，求解域内各个节点坐标、按照线性多项式确定的节点位移及节点载荷（节点等效后的荷载）罗列于下表。弹性模量取$E=1000$，泊松比$\nu =0.3$。

为简单起见，在本次模型中可取单向拉伸应力状态，位移场为：

$$\{\begin{array}{l}u=0.01x\\ v=-0.003y\end{array}$$

由以上位移场如何确定面力呢？

可以从位移场中看出$x$向应变${\epsilon }_x$设定为0.01。由泊松比$\nu =\frac{{\epsilon }_x}{{\epsilon }_y}$可推出纵向应变${\epsilon }_y=0.003$。模型的弹性模量$E=1000$，应力${\sigma }_x=E{\epsilon }_x=10$。这也就是上图所示的面力大小来源。

数值实现

试验A

对于试验A的前处理数据文件可设置为：

*Node
      1,    0.0,       0.0
      2,    2.5,       0.0
      3,    2.5,       3.0
      4,    0.0,       2.0
      5,    0.5,       0.5
      6,    2.0,       0.75
      7,    1.75,      1.75
      8,    0.65,      1.6
*Element
1, 1, 2, 6, 5
2, 2, 3, 7, 6
3, 3, 4, 8, 7
4, 4, 1, 5, 8
5, 5, 6, 7, 8
*Material
1000.0, 0.3,1.0,1
*Constr
1,1,0.
1,2,0.
2,1,0.025
2,2,0.
3,1,0.025
3,2,-0.009
4,1,0.
4,2,-0.006
5,1,0.005
5,2,-0.0015
6,1,0.02
6,2,-0.00225
7,1,0.0175
7,2,-0.00525
8,1,0.0065
8,2,-0.0048

在试验A中，所有节点的位移按照上表施加给定位移，用平面应力四节点等参单元有限元程序计算各节点力：

从上图可以看到，除了和已知的节点力一致外，其他节点力与上表对应的节点力精度保持在$10^{-15}-10^{-16}$之间，试验通过。

试验B

对于试验B的前处理数据文件可设置为：

*Node
      1,    0.0,       0.0
      2,    2.5,       0.0
      3,    2.5,       3.0
      4,    0.0,       2.0
      5,    0.5,       0.5
      6,    2.0,       0.75
      7,    1.75,      1.75
      8,    0.65,      1.6
*Element
1, 1, 2, 6, 5
2, 2, 3, 7, 6
3, 3, 4, 8, 7
4, 4, 1, 5, 8
5, 5, 6, 7, 8
*Material
1000.0, 0.3,1.0,1
*Constr
1,1,0.
1,2,0.
2,1,0.025
2,2,0.
3,1,0.025
3,2,-0.009
4,1,0.
4,2,-0.006

在试验B中，1-4节点的位移按照上表给出的位移进行施加，使用平面应力四节点等参单元有限元程序计算其余各节点位移与节点力：

从上图可以看到，各个节点力与上表对应的节点力精度保持在$10^{-15}-10^{-16}$之间，所有节点位移和节点力也是上表对应的节点位移和节点力完全吻合，试验通过。

试验C

对于试验C的前处理数据文件可设置为：

*Node
      1,    0.0,       0.0
      2,    2.5,       0.0
      3,    2.5,       3.0
      4,    0.0,       2.0
      5,    0.5,       0.5
      6,    2.0,       0.75
      7,    1.75,      1.75
      8,    0.65,      1.6
*Element
1, 1, 2, 6, 5
2, 2, 3, 7, 6
3, 3, 4, 8, 7
4, 4, 1, 5, 8
5, 5, 6, 7, 8
*Material
1000.0, 0.3,1.0,1
*Load
1,1,-10
2,1,15
3,1,10
4,1,-15
*Constr
1,1,0.
1,2,0.
2,2,0.

在试验C中，对节点1和节点2施加上图所示的边界条件在节点1、节点2、节点3和节点4的$x$方向上施加上表给定的节点载荷，使用平面应力四节点等参单元有限元程序计算各节点位移和节点力：

从上图可以看到，各个节点力与上表对应的节点力精度保持在$10^{-15}-10^{-16}$之间，所有节点位移和节点力也是上表对应的节点位移和节点力完全吻合，试验通过。

云图解释

为了更加直观的解释何为常应变状态，现将分片试验的云图绘制于下图，对于线弹性问题，若单元应力状态为常应力状态，那么该单元也为常应变状态。从下图可以看出，五个单元组成的分片区域内，应力为常量，从而表达出常应变状态。

读者也可用abaqus做abaqus自带单元的分片试验，算例文件均已给出，应力云图见下图，也是呈常应变状态，分片试验通过！

以上推文的理论内容及整套代码和算例文件已更新至《有限元基础编程百科全书》最新PDF中，存放在知识星球中，缩略图如下：

星球加入方法：

有限元分片试验 | 从理论到手搓代码step-by-step数值实现！