有限要素法、弾塑性解析についてWeb検索をしていた所、Pythonコードを見つけました。
https://qiita.com/Altaka4128/items/86e25d66196dfe7160f3
塑性域に達した時の処理について理解を深める事ができて大変ありがたいです。
実行方法
記事よりコードをコピーして以下のファイルを用意しました。
・main.py
・FEM1d.py
・Node1d.py
・d1t2.py
・Material.py
・Boundary1d.py
コマンドプロンプトでpythonコードを配置したディレクトリにいき、
py -3.12 main.py
で実行できました。
(私の環境ではpython3.9, 3.10, 3.12を共存させています。解析実行は最新のpython3.12で行いました。）
ゼロ割のWarningがでる原因はFEM1d.py内のvecDispFirstがvecDispの参照渡しとなっているからのようです。
vecDispFirst = vecDisp.copy()
としておいたほうが良さそうです。
この修正をしなくても一応解析は最後まで流れて正しい結果が得られます。これは収束判定を２つ行っているので、ゼロ割を起こす方を無視しても、もう一方の
if np.isclose(LA.norm(vecd), 0.0):
で収束判定を正しく行えているからのようです。
Abaqus Learning Edition2023で検証
記事内にAbaqusのinpファイルもあるのでAbaqusでも実行してみました。
1次元モデルなのでただの棒線ですが、変位量が一致する結果が得られました。
