广告
关于有限元应力结果精度的几点讨论
2019-02-11 09:47:39· 来源：模态空间 作者：王朋波
1、引言作为一名汽车结构CAE工程师，使用有限元法进行应力分析是我多年的日常工作。但是关于有限元应力结果的一些技术点一直未能吃透，查到的相关文献也不多，文
1、引言
作为一名汽车结构CAE工程师，使用有限元法进行应力分析是我多年的日常工作。但是关于有限元应力结果的一些技术点一直未能吃透，查到的相关文献也不多，文献内容也存在很多不尽不实之处。最近一段时间，抽空重新读了一遍王勖成老师的《有限单元法》，也查阅了一些其它资料，简单整理出了一些关于应力结果的笔记，再加上一点自己的理解，与大家分享交流。
2、有限元位移解的下限性
利用最小位能原理求得的位移近似解所对应的弹性变形能是精确解变形能的下限，即位移近似解在总体上偏小，模型偏于刚硬。利用最小余能原理得到的应力近似解所对应的弹性余能是精确解余能的上限，即应力近似解在总体上偏大，结构的计算模型偏于柔软。
我们常用的结构单元大都是位移元，以位移为未知量，基于最小位能原理建立。位移元得到的位移解具有下限性质，在给定的载荷之下，计算模型的变形比实际要小。当单元网格分割得越来越细时，位移数值解将由下方收敛于精确解，即得到真实解的下界。
需要注意的是，位移数值解并不是在每个点上都小于精确解，数值解只是总体上小于精确解，更准确的说法应该是外载荷在数值解上做的功小于在精确解上做的功，即
ũTP≤uTP （1）
上式中，ũ和u分别为位移数值解和精确解，P为外载荷。
位移解的下界性质从物理上非常好理解。连续体具有无限多个自由度，可以有无数种复杂的变形模式。划分单元后，位移场用有限个结点的位移和对应的形函数来表示，即用一些简单变形模式来逼近实际变形，如图1所示。这就意味着连续体的变形受到了约束和限制，即刚度较实际增加了。由刚度方程可知，在外力相同的情况下，所求得的位移近似解将在总体上偏小。
图1 用特定的形函数逼近准确解
3、应力结果的精度
使用位移有限元法进行结构分析时，未知的场函数是结构位移。利用最小位能原理建立的求解方程是系统的平衡方程，求解方程得到的是各结点的位移，但实际工程问题往往更关注结构应力分布。
位移有限元法求解应力的基本步骤如下：
- 引入位移边界条件；
- 求解方程组得到各结点的位移a；
- Ka=P （2）
- 根据单元各结点位移，通过导数运算求得应变和应力。
ε=Lu=Bae σ=Dε=DBae （3）
例如对于平面问题，
应变矩阵B是对插值函数求导得到的矩阵，每求导一次，插值多项式的次数就降低一次。所以通过导数运算得到的应力解的精度较位移解降低了一阶。例如线性单元的应力是近似均布的，二次单元的应力是近似线性分布的，其精度都比位移解低（虽然应力分布函数还包含了一些高次项，但这些高次项都是非完全项，并不能提高应力精度）。
有限元应力解的近似性表现在：
- 单元内部一般不满足平衡方程；
- 单元与单元交界面上应力一般不连续；
- 力的边界上一般也不满足力边界条件。
只有单元尺寸无限趋近于0时，即自由度数趋近无穷时，才能精确满足平衡方程、力边界条件和单元交界面上的应力连续性。在单元数量有限时，这些条件只能近似满足，除非实际应力变化的阶次等于或低于所用单元的应力分布函数阶次。
4、应力解的震荡性质
位移有限元法从力学上解释是求位移变分所引起的应变能为极小值的问题，从数学上解释是求解应力近似解与精确解差值的加权最小二乘问题。
与位移结果不同，位移元的应力结果并没有下限性质。应力结果是精确应力在加权最小二乘意义上的近似解，应力近似解必然在精确解上下震荡；并且在某些点上，近似解恰好等于精确解，即单元内存在最佳应力点，如图2所示。
图2 有限元应力解的震荡性质
我们在有限元强度分析中，经常发现在应力集中部位的应力解低于精确解，所以有人误以为位移有限元法的应力解也有下限性质。但真正原因是单元数目有限，且应力近似解的阶次低，无法正确描述应力的剧烈变化。对于高应力梯度区域，有限元解通常给出的是一个比实际平滑的结果。如果结构上某部位的应力突然上升，该部位的应力解将低于精确解，但如果某个部位应力突然下降，该处的应力解反而会高于精确解。
位移元的应力解没有下限性质从物理上也容易解释，有限元相当于在结构内部增加了约束，即增加了整体刚度。在外载荷相同的情况下，刚度增加能导致位移数值解偏小，但解得的应力应该在总体上保持不变，否则就无法满足平衡方程。
5、减缩积分和完全积分单元
完全积分是指当单元具有规则形状时，所用的高斯积分点的数目足以对单元刚度矩阵中的多项式进行精确积分。
减缩积分是指对单元刚度阵进行积分时，所用的高斯积分点数低于精确积分的要求，一般是按照下式来确定高斯积分点数，
n=p-m+1 （5）
式中，n为高斯积分点数，p是插值函数中完全多项式的方次，m是微分算子L中的导数阶次（对于用于弹塑性力学分析的实体单元，m=1）。
按照式（5）确定积分点数的减缩积分单元比完全积分单元在每个方向上减少了一个积分点。
对于大部分商用有限元软件，二维实体单元中的线性与二次四边形单元有减缩积分形式，积分点数分别为1*1和2*2。三维实体单元中的线性和二次六面体单元有减缩积分形式，积分点数分别为1*1*1和2*2*2。三角形单元、四面体单元和楔形单元一般没有减缩积分形式。
图3 完全积分和减缩积分单元
因为降低了刚度矩阵的积分精度，减缩积分单元看似会影响结果的准确性。但实际计算表明，采用减缩积分往往可以取得比精确积分更好的精度，主要原因如下：
1）精确积分常常是由插值函数中非完全多项式的最高方次所要求，而决定有限元精度的，通常是完全多项式的方次。这些非完全的最高方次项往往不能提高精度，反而有负面影响。采用式（5）的减缩积分方案，积分精度恰好保证完全多项式方次要求，实质是用一种新的插值函数代替原插值函数，从而改善单元精度。
2）在最小位能原理基础上建立的位移有限元，计算模型具有较实际结构偏大的整体刚度。选取减缩积分方案将使有限元模型的刚度有所降低，因此可能有助于提高精度。
3）减缩积分方案对于泛函中包含罚函数的情况也常常是必须的，以保证罚函数矩阵的奇异性。例如基于相对自由度的梁单元和板壳单元，采用完全积分会出现剪切锁死问题，导致出现完全歪曲的结果，改用减缩积分方案就可以有效解决。
4) 减缩积分因为积分点数比全积分少很多，所以大幅减少了计算量，提高了计算效率。
实际工程应用中，完全积分单元容易出现剪切锁死和体积锁死等问题，即使划分很细的网格，精度依然很差，所以一般不推荐使用。
通常应选择各种减缩积分单元或者修正单元。线性减缩积分单元只有一个中心积分点，解决了锁死问题，但存在沙漏问题，有些情况下表现的过于柔软，好在已经发展出多种有效的沙漏控制技术，只要网格比较细密，就能得到较精确的位移解。二次减缩积分单元则基本上没有沙漏与锁死问题，能够提供很高的计算精度，对于通常的变形和强度分析是比较好的选择。
6、等参单元的最佳应力点
如前述，应力近似解是应力精确解在加权最小二乘意义上的近似解。如果位移近似解是p次多项式，L是m阶微分算子，则应力近似解是p-m阶多项式。如果应力精确解是比近似解更高一阶的多项式，即p-m+1阶多项式，则在p-m+1阶高斯点上，应力近似解和精确解是相等的，即近似解在这些积分点上具有比本身高一阶的精度。这些高斯点称为单元的最佳应力点，又称优化应力点或超收敛应力点。
图4 单元最佳应力点
上面的讨论是假定了单元雅各比行列式为常数，且每个单元内的应力近似解变分独立。所以以上结论仅对结点等间距分布的一维单元是严格的，对于二维和三维单元只能是近似的。但一般情况下我们仍然能得到如下推论：
在等参元中，单元的p-m+1阶高斯点上的应力和应变比其他部位具有更高的精度，因此称p-m+1阶高斯点是单元的最佳应力点。对比式（5）可知，单元的最佳应力点恰好是减缩积分方案的积分点。
很多文献认为单元积分点就是最佳应力点，这个说法并不准确。前面关于最佳应力点的讨论并未涉及单元积分方案，无论是减缩积分还是完全积分单元，最佳应力点都是p-m+1阶高斯点。换句话说，对于减缩积分单元，最佳应力点与其积分点恰好一致，但对于完全积分单元，最佳应力点与其积分点并不一致。比如线性减缩积分四边形单元，最佳应力点就是中心积分点，但对于线性完全积分四边形单元，最佳应力点仍然是中心点，与2*2个单元积分点并不一致。
既然很多情况下，单元积分点并不是最佳应力点，为什么商用有限元软件都是计算输出积分点应力呢？主要原因是为了节省计算量和内存需求。
根据式（3）可知，计算单元内任意一点的应变或者应力，需要计算集成应变矩阵B。集成单元刚度阵时，软件已经将各积分点位置的B矩阵计算并存储，求得各结点位移后，直接调用内存中的B矩阵就能算得积分点应力。但对于非积分点位置，则需要重新计算和存储该点的B矩阵。
所以商用有限元软件为了减少计算和内存消耗，都是只计算积分点应力，然后再外推到单元其它位置。
7、减缩积分和完全积分单元的应力精度
减缩积分单元只要进行了有效的沙漏控制，就能得到较高精度的位移解，但在它的应力结果通常不能令人满意，主要表现为在应力集中区域减缩积分单元给出的应力值低于实际值。
所以有人认为减缩积分方案的位移解精度更高，而完全积分方案的应力解精度更高，其实这个观点并不正确。
既然应力结果是对位移结果求导得到的，要想得到高精度的应力解，当然要先保证高精度的位移解。减缩积分单元的位移解精度较高，而且应力输出的积分点恰好是单元最佳应力点，因此其应力精度总体上是高于完全积分单元的。
但减缩积分有一个缺点是应力输出点（即积分点）较少，无法描述单元内应力的剧烈变化，对于应力集中部位，会给出一个过分平滑的单元应力结果。尤其是线性减缩积分单元，只有一个中心应力输出点，外推的结果只能是整个单元的应力为常值。如图5所示。
图5 线性减缩积分单元的应力结果
完全积分单元在每个方向都多了一个积分点，也就是每个方向多一个应力输出点，所以在结点位移足够精确的前提下，完全积分单元能够对高应力梯度提供更好的描述。
比较好的方案应该是对结构整体上用实施沙漏控制的减缩积分单元离散，对于应力集中部位则局部采用完全积分单元。需要格外注意的是，完全积分单元只能用于不承受弯曲载荷的局部，并且要采用加密的网格，以保证位移解精度。局部采用完全积分单元的做法并不能真正提升应力解精度，但对于我们所关注的局部高应力结果有明显改善。
对于三维实体结构，高应力通常发生于结构表面，表面应力是根据单元内部积分点应力外推得到的。即使我们的结点位移是精确解，外推得到的表面应力的精度也不会很高。为克服这个缺点，可以在实体有限元模型表面覆盖一层同种材料的膜单元（建议厚度取0.01毫米）。膜单元根据表面结点位移来直接计算表面应力，避免了由实体内部积分点外推造成的误差，所以表面应力的精度将有明显改善。有些商用有限元软件已经提供了类似的表面应力计算方法，无需用户自己建立表面膜元。
8、小结
- 位移有限元法得到的位移解具有下限性质，在给定的载荷之下，计算模型的变形比实际要小。
- 位移元的应力结果并没有下限性质，应力近似解必然在精确解上下震荡。
- 等参单元的最佳应力点是p-m+1阶高斯点，与减缩积分单元的积分点一致，但与完全积分单元的积分点并不一致。
- 完全积分单元的应力精度并不会高于减缩积分单元，但在应力集中部位细化网格和局部使用完全积分单元，有助于改善局部高应力结果。
- 在实体有限元模型表面覆盖膜单元来计算表面应力，能够明显提高表面应力的精度。
作者简介
王朋波，清华大学力学博士，汽车结构CAE分析专家。重庆市科协成员、《计算机辅助工程》期刊审稿人、交通运输部项目评审专家。专业领域为整车疲劳耐久/NVH/碰撞安全性能开发与仿真计算，车体结构优化与轻量化，CAE分析流程自动化等。王朋波私人微信：poplewang。
广告
广告
最新资讯
-
NVIDIA 发布 2025 财年第三季度财务报告
2024-11-21 13:30
-
Mack卡车为买家推出创新的虚拟现场探索体验
2024-11-21 13:29
-
氢燃料电池卡车从1到100要多长时间？戴姆勒
2024-11-21 13:28
-
聚焦消费者用车极限环境，2024中国汽研汽车
2024-11-21 13:21
-
新能源汽车高寒环境可靠性行驶试验研究
2024-11-21 13:19
广告
广告
