悬臂梁解析解与有限元解对比、差异分析及“膜单元”的应用

为了对比单元类型对有限元计算结果的影响，特选取一悬臂梁做对标验证，并试图就一些现象做出解释说明。

1. 问题描述

如下图所示悬臂梁，截面为10mm*10mm的矩形，长度l=100mm，杨氏模量210GPa，泊松比0.3，受均布压力载荷为10N/mm，试求出悬臂梁端部最大扰度、约束端弯曲应力、距离约束位置l/2处的弯曲应力。（参考材料力学Ⅰ（第六版） 刘鸿文，P124-125，P149，P195）

很多资料中都会说明为什么不选择约束端应力做对标验证，大致意思就是：约束端存在应力集中，有限元计算时存在误差。本文特意提取约束端弯曲应力，看看与理论解的误差有多大。

2. 理论解：

3. 有限元解：

为了对比单元类型对有限元结果的影响，依次采用3D实体单元、2D壳单元、1D梁单元来求解结果，依次提取了自由端的最大位移、梁中部最大主应力、约束端最大主应力。

（1）一阶六面体单元（单元尺寸2mm）

（2）一阶四面体单元（单元尺寸2mm）

(3) 壳单元 （单元尺寸2mm）

(4) 梁单元（单元尺寸10mm）

在建模过程中，需要注意以下几点：

(1) 除了梁单元外，实体单元和壳单元施加均布载荷时都是通过面(face)来施加，elems>faces，因此只要输入pressure=10N/mm/10mm=1N/mm2，有限元软件中pressure单位为N/mm2或者N/m2，分母为面积单位。使用梁单元时，均布载荷是施加在1D单元上，elems>elems，pressure=10N/mm*10mm（每个单元长为10mm）=100N/单元。

(2)使用梁单元时，因为需要查看梁单元内部的应力，单元类型选用PbeamL，不能选用pbeam，否则无法查看内部应力。

4. 结果对比（假设计算误差小于5%时，认为结果正确，涂绿）

5. 结论

(1) 约束端存在严重的应力畸变，采用实体单元（不论是六面体还是四面体）时计算精度较差，但是当采用壳单元和梁单元，精度很高；

(2) 使用壳单元和梁单元，不论应变还是应力（甚至包括约束端），都能得到满足要求的结果，尽管它们的单元、节点数量很少，特别是梁单元。

6. 进一步探讨

采用3D六面体单元时，为什么应力结果相差这么大（误差在17%左右）？？

本人的愚见：在本例中，当使用解析法求解时，并没有考虑悬臂梁的厚度（厚度只是影响抗弯截面系数），能否提取3D实体单元中的面力，看看该值与解析解的差距。怎样提取3D实体单元中的面力呢？采用“膜单元”！！

在3D实体单元上构造了一层膜单元（即非常薄的壳单元），构造膜单元具体方法为：先通过face找到实体单元的表面2D，再挑选所需位置的2D单元（其余位置的删除），并赋予shell属性，材料与实体单元一样，厚度设为很小，比如1e-5mm。为了方便查看结果，把这些2D单元放到一个set中。提交求解后，查看膜单元的最大主应力，如下图。

结果显示，约束端最大主应力为299.8Mpa，中部最大主应力为72.28Mpa，与解析解已经很接近了。