有限元求解：结构应力法如何实现的网格不敏感呢？

weixin

导读：
最近看到有人分享了有限元的应力集中和应力奇异，文章写得特别精彩。在学校时我一直也很疑惑，既然应力计算结果受网格尺寸影响，那到底有限元是否有能力识别应力集中呢？

这个问题直到接触到董平沙教授提出的“结构应力法”才得以解决。结构应力是董平沙教授研究焊缝疲劳开裂机理时，针对焊接结构疲劳问题的特殊性并基于力学原理而定义的一种应力，有兴趣的可以查阅相关文献。写篇介绍性的文章，让大家了解下结构应力方法。

01、悬臂梁的计算困惑

下面是一个受弯悬臂梁的例子，如图1所示，该梁尺寸为10mm×2mm×1mm，根据材料力学公式，梁的中部上下表面两个点的理论解为±750MPa。

图1 悬臂梁示意图

图2 悬臂梁的两单元计算结果

根据几何模型建立有限元模型，单元类型为平面应力的PLANE182单元，划分两个单元，求解后，单独选取单元2，则中部节点3和6的应力分别为±117.6MPa，这个结果与理论解相差很远，网格细化后，应力肯定会逐渐升高，感兴趣的小伙伴可以试一试到底细化到什么程度能和理论解对上。

这个例子本身没有什么意义，但却说明一个很重要的问题，有限元求解的值是关于网格敏感的。对于一个新手而言，不知道所在行业的结构划分多大尺寸的网格才能和试验结果相对应，这是很可怕的。

02、如何实现的网格不敏感

董平沙教授提出的结构应力法如何实现的网格不敏感呢？想要知道具体理论的小伙伴可以参考书籍《焊接结构抗疲劳设计理论与方法》。咱们还是说上面的简支梁（这个例子也是董老师的，借来用用）。

在有限元计算时，单元边上的分布载荷要向节点转化，这个学过有限元分析都应该了解，在有限元原理书籍中也叫等效载荷。而董老师网格不敏感的结构应力法利用节点力求线力和线矩时，却是这个过程的逆过程，即有限元求得的节点力和力矩转化为线力和线矩，该线力与线矩是指焊线处单位长度上的力与力矩。具体怎么转化的呢？

图3 悬臂梁的两单元计算结果

如图3所示，节点1和2在y 轴方向的节点力及绕x 轴的力矩分别为Fy1、Fy2 和Mx1、Mx2；y 轴方向单元边的线力及绕x 轴的线力矩分别为fy1，fy2 和mx1、mx2 根据力的平衡方程，可以求得下式（小伙伴可以用材料力学公式对节点1列平衡方程，即求合力及合弯矩）。

求上等式右侧项矩阵的逆可得：

所以节点1及节点2的结构应力定义式：

再来看那个简支梁，将单元2隔离出来后，节点3和节点6的单元节点力Fy 分别为±250N，没有节点弯矩，如图3所示，在图3中还用箭头描述了弯曲情况下截面的应力分布。节点3和节点6的距离为l=2mm，板厚t=1mm，根据上式的节点3和节点6的结构应力为±750MPa，这个结构应力与理论结果一致。

图3 悬臂梁的两单元节点力计算结果

通过该方法给出了与梁理论完全一致的结果，这是在意料之中的，因为梁是在静定载荷条件之下。该方法的中平衡理论，要求这样计算出的应力应当满足平衡条件。需要注意的是，这个例子仅仅用来阐述该方法的基本概念，而不是表示该方法可以没有任何限制的使用，设计的建模和网格问题，有兴趣的小伙伴可以查阅相关资料。

基于有限元的基本理论，网格不敏感的结构应力法的力学解释并不困难。结构应力的计算基础是认为截面上节点力的合力一定与外力平衡，这样在划分有限元网格时，假如在给定外力的前提下划分20个网格的节点力的合力与这个外力平衡，那么划分10个节点力的合力也应与这个外力平衡。

事实上，网格不敏感的原因完全可以从力的平衡角度去解释，因为在外力被给定的情况下，同样一条焊趾或焊线上的节点力的个数无论多少，其合力都将与这个外力平衡，所以结构应力对有限元网格一定不敏感。当创建有限元模型时，需要注意的仅仅是在基于虚功原理将节点力转换为线力或线力矩时要有足够的节点即可。

注意，结构应力的网格不敏感与通常的应力计算过程中为了提高计算精度而加密网格是完全不同的两个思考方向，由于在计算结构应力时追求过细的网格是没有必要的，所以结构应力对有限元网格不敏感这一特点在应用过程中有重要价值，因为这将显著降低有限元网格数量，从而降低了对计算机硬件能力的要求，在硬件配置不是很高的计算机上，也可以完成较大规模的计算任务。

03、实际结构中的应用

下面看看具体怎么应用到焊缝中的，图4给出了外力作用下在焊缝截面上沿着厚度方向的应力分布，这个应力分布因含有缺口应力而呈现出高度非线性。

虽然理论上或数值上直接求解这个非线性应力分布是很困难的，但是我们可以将这个高度非线性的应力进行分解。分解后，第一部分应力是只与外力相关，且与外力互相平衡的那一部分；第二部分应力则是去掉了第一部分而余下的应力分布。它虽然包含了非线性部分，但是由于第一部分已经与外力平衡，因此这部分称之为缺口应力的分布一定处于自平衡状态

图4 焊趾截面内的应力分解示意图

有了上述分析就可以定义第一部分为结构应力，第二部分为缺口应力，二者之和，就是原来截面上非线性应力的分布。这里忽略了剪切应力的影响，一个原因是它对焊缝开裂的贡献很小，如果剪切应力不能被忽略，那将是多轴疲劳问题，这里暂不讨论。那么，为什么第一部分应力称之为结构应力呢？因为，它是满足平衡条件并可以用结构力学的方法计算得到。

不失一般性，将一段焊缝划分成n 个单元时，节点编号由1至n，焊线上的各节点距离为l1 至ln-1，根据力的平衡方程，可求得各节点力Fyn 与线力fyn 的对应关系。

上式中矩阵L 只与节点距离相关，这里定义为单元长度等效矩阵：

采用Matlab等计算软件可对L 矩阵求逆，结果为L-1，由于求逆后数据不再是带状，这里不详细列出，感兴趣的小伙伴试一试。

由此可得线力与节点力的对应关系：

同理，线矩mx与节点力矩Mx 与上述表达式形式相同，这样当有n 个节点在相同的单元厚度t 的情况下，各节点的结构应力σn，以及各节点力Fyn 及力矩Mxn 可用矩阵方式表为：

04、算例分析

迄今为止，已经通过董老师结构应力方法计算了很多算例，这里也不过多介绍，引用董平沙老师比较经典的案例给大家展示一下这种方法的优越性。下面是三种单元尺寸不同的接头，长度分别为0.5t、1t 和2t，分别对比最大主应力、Von Mises 应力和结构应力，只有结构应力结构基本保持不变，最大主应力和Von Mises 应力随着网格尺寸的增加而增大。

图5 三种不同有限元网格尺寸计算的三类应力结果对比

必须声明，文章并非原创，都是引用了董平沙教授、兆文忠教授、李向伟博士等的研究成果，因为感觉很多人不知道这种方法，所以单独将结构应力计算方法进行了整理出来，通过简单的内容介绍给大家，希望有兴趣的参阅《焊接结构抗疲劳设计理论与方法》。

来源：仿真秀App微信公众号（ID：fangzhenxiu2018），作者：云兵。文章略有删减。