有限元编程深度解析：节点应力计算与磨平处理技巧
在应用位移有限元法分析时，直接的未知量是节点位移，实际问题中，可能最为关心的是结构体的应力分布，尤其是最大应力的位置和数值。需要根据以下公式由内核求解程序计算得到的节点位移推算出单元内的应力。
应变关系矩阵 ，是形函数矩阵对坐标进行求导后得到的矩阵。求导一次，形函数的阶次就降低一次，进而通过求导运算得到的应变、应力就较节点位移精度降低一些。
此时，延伸出一个概念：最佳应力点。在张雄教授的《有限元法基础》和王勖成教授的《有限单元法》中均有提及，感兴趣的小伙伴也可以在谷歌学术中搜索相关的文献进行学习。
本期将围绕应力问题展开分析，从节点应力的求解到共节点应力磨平处理，让木木带你回味这些曾经枯燥的知识点。
主要包括以下内容：
节点应力计算方法 高斯积分点外插法 直接坐标法 应力磨平 直接绕节点平均 面积/体积加权平均 Abaqus 内部应力磨平方法（avg=75%）
节点应力计算
单元内高斯积分点处的应力精度较高，但是在后处理时常常用到的是节点处的应力，同时节点处的应力也是更加直观，实用性更强。一般情况下，节点处的应力大于结构内部。本小节主要介绍两种常见处理节点应力的方法—高斯积分点外插、直接坐标法。
高斯积分点外插
以二维等参元为例，高斯积分点处的应力记为 ，节点处的应力记为 ， 表示单元形函数在第一个高斯积分点处的第二个形函数值，有：
对上式求逆，可得：
其中， ， ，
同理，对于三维八节点等参元，有：
其中， ， ， ， 。
以上的过程就是传说中的“应力外插法”。如果要开发通用性较强的程序，需要预先准备好各种单元类型的型函数逆矩阵，用点心，其实也不是很麻烦。对于像我这种懒人，自然想再看看有没有更方便的处理方法，下面一节将会给大家带来非常简单粗暴的方法！
直接坐标法
用户也可以直接将单元节点处的坐标值代入推文开头的应力计算公式，也同样可以得到节点处的应力，不过如果节点坐标中包含 0 元素，则应力应变的值会出现"NAN"，在之前的推文验证轴对称单元时出现此类问题，做了一点人为的操作，详细请看：直接坐标法修正。
我在程序中是这样处理的，仅供参考：
[x, e] = obj.localCoord2D(); % 获取单元节点坐标值
elenode = length(elementNode);
for n = 1:elenode
[B, Jacob] = Bmatrix2D(elenode, elemNodeCoordinate, obj.DOF, x(n), e(n)); % 将单元节点坐标直接带入公式进行计算
sigma = D * B * delta;
...
end
应力磨平
经过上面两种方式，都可以获取每个单元的各个节点处的应力值，但是问题就随之而来了，对于共节点处的应力该怎么处理呢？请看下面的示意图：
对于中间共节点处，共有四个单元，每个单元在该节点处的应力都不相同，那么共节点处的应力是到底该怎么表达呢？如何由左图的应力差异性转变为右图的应力平滑呢？
直接绕节点平均
比较简单粗暴的办法就是把该节点处的所有单元节点应力累加起来除以单元数，即： ，规范的数学表达：
其中， 表示共节点 上的平均应力， 表示各个单元计算所得到的在共节点 上的应力， 为围绕该共节点周围的全部单元。
面积、体积加权平均
当单元面积/体积相差不大时，“直接绕节点平均”可以得到比较合理的结果，但是相邻单元如果面积/体积相差比较大时，应力的近似解误差会比较大，这时可以将面积/体积考虑在内，对其进行加权平均，即：
规范的数学表达：
注：以上处理只是计算后处理中的一种局部改善，并不能从根本上解决节点应力、应变精度差的问题。
这里给出已有的代码，来自《结构分析的有限元法与 MATLAB 程序设计》，仅供参考。
gNodeStress = zeros( node_number, 6 ) ;
for i=1:node_number
disp( sprintf( '计算节点应力，当前结点: %d', i ) ) ;
S = zeros( 1, 3 ) ;
A = 0 ;
for ie=1:1:element_number
for k=1:1:3
if i == gElement( ie, k )
area= ElementArea( ie ) ;
S = S + gElementStress(ie,1:3 ) * area ;
A = A + area ;
break ;
end
end
end
gNodeStress(i,1:3) = S / A ;
gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
end
Abaqus 内部应力磨平方法
相信习惯用 Abaqus 用户对 并不感到陌生吧，那它到底什么含义呢？在程序中又该如何实现呢？
这里推荐吉林大学左文杰教授在 B 站上的有限元编程公开课，详细讲述了应力磨平时阈值的设置原理，视频讲的很好，我都觉得没必要再写一遍了，大家可直接看原视频的讲解。
以上就是给大家带来的有限元编程时节点应力的处理方法，希望可以帮助到有需要的小伙伴，感谢您的阅读，我们下期再见~
参考文献
[1]. Cook, Robert D. Concepts and applications of finite element analysis. John wiley & sons, 2007.
[2]. 王勖成. 有限单元法. 清华大学出版社有限公司, 2003.
[3]. 曾攀. 有限元分析基础教程. 北京: 清华大学出版社, 2008.
[4]. 徐荣桥. 结构分析的有限元法与 MATLAB 程序设计. 人民交通出版社, 2006.
