上一篇技文《ABAQUS结果提取大于某值的区域体积-CAE方法》中带大家使用ABAQUS CAE界面直接提取大于100e6 Mises等效应力的区域体积,方法虽然比较好操作,但也存在明显的缺点:
方法太繁琐;统计历程曲线时会显得非常麻烦,因此我们找到了代步工具:Python。
精度较低;这是由于一个单元内只需要一个积分点满足数值要求,整个单元都会显示出来,但一个单元可能存在多个积分点,比如演示案例中采用六面体二阶减缩积分单元(C3D20R)存在8个积分点,单元内部分积分点可能不满足条件,也会被统计在内,从而导致统计的体积是偏大,我们对比下CAE方法和Python方法的结果差异如下图,最少都差32%!
样的误差虽然可以通过提高网格密度来减缓,但并不能完全避免,因此也就有必要对这些单元内所有积分点值进行判断,获得单元内满足条件的积分点所占比例,再对其体积进行加权相加。而这样的操作过程显然不适合手动统计,再一次将目光投向了我们的代步工具:Python!
没有Python基础的小伙伴,建议先看曹金凤姐姐的《Python语言在Abaqus中的应用》或江丙云哥哥的《ABAQUS Python二次开发攻略》,不然会有些吃力。虽然代码备注的已经非常详细,但还是需要一些背景支撑的。
编程需求与分解
为了与CAE过程进行对比,我们还是将需求定义为:提取悬臂梁加载过程中Mises等效应力大于100e6的区域体积。
悬臂梁模型背景:注意采用的是C3D20R单元,每个单元拥有20个节点、8个积分点;(不懂积分点含义的需要补充些有限元理论了)
一级需求分解:
1.遍历所有Step分析步和Frame帧内的Mises场输出,判断是否大于100e6;
2.将满足条件的单元体积进行累加;
3.求解每个单元的体积;
二级需求分解:
1.1遍历Step
1.2遍历Frame
1.3遍历Mises场输出
1.4判断数值是否大于100e6
2.1需要考虑:部分单元内积分点n<=8个满足判断条件的情况。可将单元体积均分为8份,当判断单个积分点满足要求时,增加1/8个单元体积;
3.1 C3D20R单元为六面体,求解任意凸六面体单元体积:可将六面体切为以6个面为地面,内部任意一点(中心点)为顶点的6个四棱锥;每个四棱锥可分为2个四面体,分别求四面体体积再相交,则为整个六面体体积;(由于是二阶单元,如果需要更高精度,可以将单个四棱锥刨分为更多的四面体)
3.2 获得六面体内部任一点,(真理:凸六面体八个角点坐标的平均值(中心点)一定位于该六面体内部),如果不是凸六面体,单元刚度将为负,所以不存在;
3.3 求任意四面体体积;
3.4 获得六面体的顶点坐标:C3D20R的前8个节点编号为六面体角点编号。
部分代码
通过上面的二次需求分解,基本将一个看似复杂的需求划分为若干个小而简单的需求,而每个小的需求通常使用十行以内代码即可完成;大家可以动笔自己编写一下,这里就不展示最终代码了。