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

一、导读

本人从大一时，上第一节力学课就开始接触有限元，那个时候老师曾说：“学会我讲的这门有限元，我可以保证你在毕业时候一定找到工作”。那个时候有限元很火，很多高校都开了课，这一晃都过去13年啦。最近看到有人分享了有限元的应力集中和应力奇异，文章写得特别精彩。在学校的时候我一直也很疑惑，既然应力的计算结果受网格尺寸影响，那到底有限元能不能有能力识别应力集中呢？这个问题直到接触到董平沙教授提出的“结构应力法”才得以解决。结构应力是董平沙教授研究焊缝疲劳开裂机理时，针对焊接结构疲劳问题的特殊性并基于力学原理而定义的一种应力，有兴趣的可以查阅相关文献。写篇介绍性的文章，让大家了解下结构应力方法，欢迎各位亲对该方法进行深入的了解，扩展其应用范围，提出不同的看法或讨论，欢迎来微信群交流。

二、悬臂梁的计算困惑下面是一个受弯悬臂梁的例子，如图1所示，该梁尺寸为

，根据材料力学公式，梁的中部上下表面两个点的理论解为

。

根据几何模型建立有限元模型，单元类型为平面应力的PLANE182单元，划分两个单元，求解后，单独选取单元2，则中部节点3和6的应力分别为

，这个结果与理论解相差很远，网格细化后，应力肯定会逐渐升高，感兴趣的小伙伴可以试一试到底细化到什么程度能和理论解对上。

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

三、 如何实现的网格不敏感董平沙教授提出的结构应力法如何实现的网格不敏感呢？想要知道具体理论的小伙伴可以参考书籍《焊接结构抗疲劳设计理论与方法》。咱们还是说上面的简支梁（这个例子也是董老师的，借来用用）。在有限元计算时，单元边上的分布载荷要向节点转化，这个学过有限元分析都应该了解，在有限元原理书籍中也叫等效载荷。而董老师网格不敏感的结构应力法利用节点力求线力和线矩时，却是这个过程的逆过程，即有限元求得的节点力和力矩转化为线力和线矩，该线力与线矩是指焊线处单位长度上的力与力矩。具体怎么转化的呢？

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

求（1）等式右侧项矩阵的逆可得：

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

再来看那个简支梁，将单元2隔离出来后，节点3和节点6的单元节点力

分别为

，没有节点弯矩，如图3所示，在图3中还用箭头描述了弯曲情况下截面的应力分布。节点3和节点6的距离为

，板厚

，根据公式（3）的节点3和节点6的结构应力为

，这个结构应力与理论结果一致。

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

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

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

四、实际结构中的应用下面看看具体怎么应用到焊缝中的，图4给出了外力作用下在焊缝截面上沿着厚度方向的应力分布，这个应力分布因含有缺口应力而呈现出高度非线性。虽然理论上或数值上直接求解这个非线性应力分布是很困难的，但是我们可以将这个高度非线性的应力进行分解。分解后，第一部分应力是只与外力相关，且与外力互相平衡的那一部分，第二部分应力则是去掉了第一部分而余下的应力分布。它虽然包含了非线性部分，但是由于第一部分已经与外力平衡，因此这部分称之为缺口应力的分布一定处于自平衡状态。

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

至

，根据力的平衡方程，可求得各节点力

与线力

的对应关系。

（4）上式中矩阵

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

采用Matlab等计算软件可对

矩阵求逆，结果为

，由于求逆后数据不再是带状，这里不详细列出，感兴趣的小伙伴试一试。

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

同理，线矩mx与节点力矩Mx与上述表达式形式相同，这样当有n个节点在相同的单元厚度t的情况下，各节点的结构应力

，以及各节点力

及力矩

可用矩阵方式表为：

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

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

感兴趣的朋友，可以加入仿真秀力学与有限元分析交流群，预计在11-12月，我们将围绕网格不敏感这个主题，将邀请本文作者云兵老师和神秘嘉宾在我们的微信群开展超级话题讨论，欢迎大家来一起讨论和学习。

作者：云兵老师，仿真秀专栏

作者声明：原创文章，首发仿真秀，部分图片源自网络，如有不当请联系我们，欢迎分享，禁止私自转载，转载请联系我们。