[ASA-4] 使用Python脚本读取单元应力并绘制莫尔圆

前不久，看博士同学在朋友圈分享李曙光老师的《复合材料损伤、破坏领域理性安在？》。其中有一页PPT戏谑研究复合材料Mohr强度准则，并点评“Mohr准则并没有那么妇孺皆知！”

Mohr准则在复合材料领域是否妇孺皆知，我不清楚，恕我孤陋寡闻。莫尔圆在金属应力分析中用几何法计算主应力，我还是略懂的。下面讲一下使用几何法画莫尔圆，计算主应力。当然，“明知故算”，帮助理解，为此写了两端Python脚本，修改参数，即画莫尔圆。

先回顾一下主应力和莫尔圆。注意下面两图中的φ和2φ。

建个简单模型，进行静力分析。材料：铝合金。孔边单元集合Em。

看一下分析结果中几种应力。对于各向同性材料，没有指定材料方向，因此，S11、S22、S12为默认全局坐标系下三个应力分量。第二行3个图从左至右依次为最大面内主应力、最小面内主应力、Tresca应力（第三强度理论，最大主应力与最小主应力的差）。 

使用下面Python脚本把圆孔周边单元的应力写出来看一下。下图中，maxIPP，面内最大主应力；maxP，最大主应力。 

上图数据文件名Data-S.txt，从ODB文件读取数据写入Data-S.txt文件的Python脚本如下。

使用几何法创建莫尔圆的Python脚本如下。下图中，n1、n2对应Data-S.txt文件中第n1行至n2-1行，每行为1个单元的应力数据。根据每行应力数据为对应单元绘制一个莫尔圆Part。Part名称后缀数字为单元编号。

查看其中一个莫尔圆，左右两端的竖线使用ABAQUS计算的面内最大、最小主应力的值绘制，即maxIPP、minIPP。圆使用S11、S22、S12三个应力分量绘制。如果竖线与圆相切，说明莫尔圆绘制正确。下图为45号单元的莫尔圆。

测一下角度，2φ=24.11°，见下图。

查看P1、P2点坐标见下图，P1→P2方向即为该单元最大面内主应力方向，P1P2与X轴夹角为φ。

显示单元45的应力S11，如下图所示。

以P1点为原点，P2点为X轴创建一个直角坐标系CSYS-1，如下图左图。将单元45的应力S11转换为CSYS-1坐标系下，如下图右上所示。显示单元45的最大面内主应力，如下图右下所示。下图中右侧上下两图应力基本一致，差异小于0.2%，说明坐标转换后S11在主应力方向上，即CSYS-1的X轴向为主应力方向。

莫尔圆可能不是妇孺皆知。但经过以上操作，应该清楚莫尔圆是咋回事了。