星辰技文|ABAQUS结果提取大于某值的区域体积-Python方法

上一篇技文《ABAQUS结果提取大于某值的区域体积-CAE方法》中带大家使用ABAQUS CAE界面直接提取大于100e6 Mises等效应力的区域体积，方法虽然比较好操作，但也存在明显的缺点：

1. 方法太繁琐；统计历程曲线时会显得非常麻烦，因此我们找到了代步工具：Python。

2. 精度较低；这是由于一个单元内只需要一个积分点满足数值要求，整个单元都会显示出来，但一个单元可能存在多个积分点，比如演示案例中采用六面体二阶减缩积分单元（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个节点编号为六面体角点编号。

 通过上面的二次需求分解，基本将一个看似复杂的需求划分为若干个小而简单的需求，而每个小的需求通常使用十行以内代码即可完成；大家可以动笔自己编写一下，这里就不展示最终代码了，需要下载的小伙伴可以点付费下载，收取少量稿费，还望大家谅解与支持。

【收费下载文件】

Lug_614.cae：为测试的模型文件，高于ABAQUS6.14软件版本则可直接打开；

Lug.inp：或可直接导入inp文件进行计算，也可abaqus command下提交运算

Lug.odb：由Lug_614.cae的job或Lug.inp计算获得，运行代码前请先计算模型。

20181225GetVolumeOfMisesGreaterThan.py：为主题相关的python代码

【义务与合作声明】

1、购买代码原则上不提供答疑，如需修改程序或定制程序，可以通过技术邻答疑和项目方式寻求合作；

2、由于工作原因，答疑预约时间为休息日和工作日的晚上8点到10点，项目合作仅能预约周末节假日；

3、预约答疑的小伙伴可在技术邻直接发送答疑和项目申请给【星辰_北极星】，看到后会立马给予回复，超出本人能力范围的将直接退还预付款。