概要
この記事では、
・Total Lagrangian法
・Green-Lagrangeひずみ
・第2Piola-Kirchhoff応力
・Saint-Venant Kirchhoffモデル
を用いた3次元有限変形FEMソルバーをPythonで実装します。
最終的には片持ち梁の大変形解析を実行し、 理論解との比較を行います。
前置き
前回の記事では微小変形を前提とした弾塑性解析プログラムを実装してみました。しかし、実際の解析では、弾塑性を考慮するなら大変形を扱いたいケースも多く存在します。そこで、大変形を考慮した解析プログラムを実装してみたいと思います。ただ、いきなり大変形+弾塑性を実装して理解するのは難しいので、まずは大変形の弾性解析を実施して、仕組みを理解してみます。
有限変形理論のための基礎知識
変形勾配テンソル
変形勾配テンソルは下記の式で定義されるテンソルです。
\begin{align}
&\boldsymbol{F} = \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{X}}
=
\begin{bmatrix}
\frac{\partial x_1}{\partial X_1} && \frac{\partial x_1}{\partial X_2} && \frac{\partial x_1}{\partial X_3} \\
\frac{\partial x_2}{\partial X_1} && \frac{\partial x_2}{\partial X_2} && \frac{\partial x_2}{\partial X_3} \\
\frac{\partial x_3}{\partial X_1} && \frac{\partial x_3}{\partial X_2} && \frac{\partial x_3}{\partial X_3}
\end{bmatrix}
\\
&X_i : 変形前の座標の成分 \\
&x_i : 変形後の座標の成分
\end{align}
ここで、$\boldsymbol{X}$は変形前の物体内の座標、$\boldsymbol{x}$は変形後の物体内の座標で、下記の図のように表されます。
変形勾配テンソルは下記のように分解することができます。
\begin{align}
&\boldsymbol{F} = \boldsymbol{V}\boldsymbol{R} = \boldsymbol{R}\boldsymbol{U} \\
&\boldsymbol{R} : 回転行列 \\
&\boldsymbol{V} : 左ストレッチテンソル \\
&\boldsymbol{U} : 右ストレッチテンソル \\
\end{align}
この分解は極分解と呼ばれます。
この変形勾配テンソルを使ってひずみを表します。
グリーン・ラグランジュひずみ
微小ひずみは物体が変形せずに回転する場合、変形していないにも関わらず、ひずみの値を持ってしまいます。
そこで、新しく物体が回転しても影響が無いひずみとしてグリーン・ラグランジュひずみを定義します。
グリーン・ラグランジュひずみ$\boldsymbol{E}$は下記の式で定義されます。
\boldsymbol{E} = \frac{1}{2}(\boldsymbol{F}^T\boldsymbol{F} - \boldsymbol{I})
ここで、$\boldsymbol{F}^T\boldsymbol{F}$は極分解により下記のように表せるため、回転成分を打ち消すことができます。
\boldsymbol{F}^T\boldsymbol{F} = (\boldsymbol{R}\boldsymbol{U})^T \boldsymbol{R}\boldsymbol{U} = \boldsymbol{U}^T \boldsymbol{U}
ここで、グリーン・ラグランジュひずみの各成分は下記のように表せます。
なお、記載にはアインシュタインの総和規約を用いています。
\begin{align}
&\boldsymbol{E}_{ij} = \frac{1}{2}\{\frac{\partial u_i}{\partial X_j} + \frac{\partial u_j}{\partial X_i} + (\frac{\partial u_k}{\partial X_i} \frac{\partial u_k}{\partial X_j})\} \\
&u_i : 変位の成分
\end{align}
第2ピオラ・キルヒホッフ応力
コーシー応力(真応力)を$\boldsymbol{\sigma}$とすると第2ピオラ・キルヒホッフ応力$\boldsymbol{S}$は下記のように定義されます。
\begin{align}
&\boldsymbol{S} = J \boldsymbol{F}^{-1} \boldsymbol{\sigma} \boldsymbol{F}^{-T} \\
&J = \det \boldsymbol{F}
\end{align}
第2ピオラ・キルヒホッフ応力はグリーンひずみと対になる(仕事共役)応力です。
仕事共役とは、仮想仕事の原理で使われる仮想仕事$\delta W$を下記のように表すことができる応力とひずみの組み合わせのことです。ここで、$:$はテンソル同士の内積を表しています。
\delta W = \int_V \boldsymbol{S} : \delta \boldsymbol{E} dV
トータルラグランジュ法
仮想仕事の原理
仮想仕事の原理により第２ピオラ・キルヒホッフ応力とグリーンひずみの関係は、下記のように表されます。
\begin{align}
&\int_V \boldsymbol{S}:\delta \boldsymbol{E}dV = \delta \boldsymbol{R}
\tag{1}
\\
&\delta \boldsymbol{E} : 仮想グリーン・ラグランジュひずみ \\
&\delta \boldsymbol{R} : 外力の仮想仕事
\end{align}
式(1)の左辺を成分毎に表すと、下記のようにベクトル形式で置き換えることが可能です。
\begin{align}
&\int_V \boldsymbol{S}:\delta \boldsymbol{E}dV = \int_V S_{ij}\delta E_{ij} dV = \int_V \{\delta \boldsymbol{E}\}^T \{\boldsymbol{S}\} dV
\tag{2}
\\
&\{\boldsymbol{S}\} =
\begin{bmatrix}
S_{11} \\
S_{22} \\
S_{33} \\
S_{12} \\
S_{23} \\
S_{13}
\end{bmatrix}
,
\{ \delta \boldsymbol{E}\} =
\begin{bmatrix}
\delta E_{11} \\
\delta E_{22} \\
\delta E_{33} \\
2 \delta E_{12} \\
2 \delta E_{23} \\
2 \delta E_{13}
\end{bmatrix}
\end{align}
$\delta \boldsymbol{E}$の成分はグリーンひずみの定義式より、下記のように表すことができます。
\begin{align}
\delta E_{ij} = \frac{1}{2}\{\frac{\partial \delta u_i}{\partial X_j} + \frac{\partial \delta u_j}{\partial X_i} + (\frac{\partial \delta u_k}{\partial X_i} \frac{\partial u_k}{\partial X_j})\ + (\frac{\partial u_k}{\partial X_i} \frac{\partial \delta u_k}{\partial X_j})\} \tag{3}
\end{align}
6面体要素内の変位は形状関数により下記のように表すことができます。
\begin{align}
& u_i = N_n u_{n,i} \\
\tag{4}
\\
&N_n : 形状関数 \\
&u_{n,i} : n節点の変位の成分i
\end{align}
$\delta u_i = N_n \delta u_{n,i}$より、式(4)は下記のようになります。
\begin{align}
\delta E_{ij} =& \frac{1}{2}\{ \frac{\partial N_n}{\partial X_j} \delta u_{n,i} + \frac{\partial N_n}{\partial X_i} \delta u_{n,j} +
\frac{\partial N_n}{\partial X_i} \frac{\partial u_k}{\partial X_j} \delta u_{n,k} + \frac{\partial u_k}{\partial X_i} \frac{\partial N_n}{\partial X_j} \delta u_{n,k}\}
\tag{5}
\end{align}
ここで、$\delta \boldsymbol{E}$をベクトル形式にすると下記のように表すことができます。
\begin{align}
&\{ \delta \boldsymbol{E}\} = \boldsymbol{B}_{L0} \{ \delta \boldsymbol{u}\} + \boldsymbol{B}_{L1} \{\delta \boldsymbol{u}\}
\tag{6}
\\
&\boldsymbol{B}_{L0} =
\begin{bmatrix}
\frac{\partial N_1}{\partial X_1} && 0 && 0 && \dots && \frac{\partial N_8}{\partial X_1} && 0 && 0 \\
0 && \frac{\partial N_1}{\partial X_2} && 0 && \dots && 0 && \frac{\partial N_8}{\partial X_2} && 0 \\
0 && 0 && \frac{\partial N_1}{\partial X_3} && \dots && 0 && 0 && \frac{\partial N_8}{\partial X_3} \\
\frac{\partial N_1}{\partial X_2} && \frac{\partial N_1}{\partial X_1} && 0 && \dots && \frac{\partial N_8}{\partial X_2} && \frac{\partial N_8}{\partial X_1} && 0 \\
0 && \frac{\partial N_1}{\partial X_3} && \frac{\partial N_1}{\partial X_2} && \dots && 0 && \frac{\partial N_8}{\partial X_3} && \frac{\partial N_8}{\partial X_2} \\
\frac{\partial N_1}{\partial X_3} && 0 && \frac{\partial N_1}{\partial X_1} && \dots && \frac{\partial N_8}{\partial X_3} && 0 && \frac{\partial N_8}{\partial X_1}
\end{bmatrix}
\\
&\boldsymbol{B}_{L1} =
\begin{bmatrix}
\frac{\partial u_1}{\partial X_1}\frac{\partial N_1}{\partial X_1} && \frac{\partial u_2}{\partial X_1}\frac{\partial N_1}{\partial X_1} && \frac{\partial u_3}{\partial X_1}\frac{\partial N_1}{\partial X_1} && \dots && \frac{\partial u_3}{\partial X_1}\frac{\partial N_8}{\partial X_1} \\
\frac{\partial u_1}{\partial X_2}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_2}{\partial X_2}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_3}{\partial X_2}\frac{\partial N_1}{\partial X_2} && \dots && \frac{\partial u_3}{\partial X_2}\frac{\partial N_8}{\partial X_2} \\
\frac{\partial u_1}{\partial X_3}\frac{\partial N_1}{\partial X_3} && \frac{\partial u_2}{\partial X_3}\frac{\partial N_1}{\partial X_3} && \frac{\partial u_3}{\partial X_3}\frac{\partial N_1}{\partial X_3} && \dots && \frac{\partial u_3}{\partial X_3}\frac{\partial N_8}{\partial X_3} \\
\frac{\partial u_1}{\partial X_2}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_1}{\partial X_1}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_2}{\partial X_2}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_2}{\partial X_1}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_3}{\partial X_2}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_3}{\partial X_1}\frac{\partial N_1}{\partial X_2} && \cdots && \frac{\partial u_3}{\partial X_2}\frac{\partial N_8}{\partial X_1} + \frac{\partial u_3}{\partial X_1}\frac{\partial N_8}{\partial X_2} \\
\frac{\partial u_1}{\partial X_2}\frac{\partial N_1}{\partial X_3} + \frac{\partial u_1}{\partial X_3}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_2}{\partial X_2}\frac{\partial N_1}{\partial X_3} + \frac{\partial u_2}{\partial X_3}\frac{\partial N_1}{\partial X_2} && \frac{\partial u_3}{\partial X_2}\frac{\partial N_1}{\partial X_3} + \frac{\partial u_3}{\partial X_3}\frac{\partial N_1}{\partial X_2} && \dots && \frac{\partial u_3}{\partial X_2}\frac{\partial N_8}{\partial X_3} + \frac{\partial u_3}{\partial X_3}\frac{\partial N_8}{\partial X_2} \\
\frac{\partial u_1}{\partial X_3}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_1}{\partial X_1}\frac{\partial N_1}{\partial X_3} && \frac{\partial u_2}{\partial X_3}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_2}{\partial X_1}\frac{\partial N_1}{\partial X_3} && \frac{\partial u_3}{\partial X_3}\frac{\partial N_1}{\partial X_1} + \frac{\partial u_3}{\partial X_1}\frac{\partial N_1}{\partial X_3} && \cdots && \frac{\partial u_3}{\partial X_3}\frac{\partial N_8}{\partial X_1} + \frac{\partial u_3}{\partial X_1}\frac{\partial N_8}{\partial X_3} \\
\end{bmatrix}
\\
&\{ \delta \boldsymbol{u}\} =
\begin{bmatrix}
u_{1,1} \\
u_{1,2} \\
u_{1,3} \\
\vdots \\
u_{8,3} \\
\end{bmatrix}
\end{align}
ここで、以下のように$\boldsymbol{B}_{L}$を定義すると、
\boldsymbol{B}_{L} = \boldsymbol{B}_{L0} + \boldsymbol{B}_{L1} \tag{7}
式(6)は以下のように表すことができます。
\{ \delta \boldsymbol{E}\} = \boldsymbol{B}_{L} \{ \delta \boldsymbol{u}\} \tag{8}
すると、式(2)は以下のように表すことができます。
\int_V \boldsymbol{S}:\delta \boldsymbol{E}dV = \int_V \{\delta \boldsymbol{E}\}^T \{\boldsymbol{S}\}dV = \{ \delta \boldsymbol{u}\}^T \int_V \boldsymbol{B}_{L}^T \{\boldsymbol{S}\} dV = \delta \boldsymbol{R} \tag{9}
ここで、式(9)の右辺は下記のように置くことができるので、
\begin{align}
&\delta \boldsymbol{R} = \{ \delta \boldsymbol{u}\}^T \{\boldsymbol{R}\} \tag{10} \\
&\{\boldsymbol{R}\} : 外力ベクトル
\end{align}
下記の式が成り立ちます。
\{\boldsymbol{R}\} = \{\boldsymbol{Q}\} = \int_V \boldsymbol{B}_{L}^T \{\boldsymbol{S}\} dV \tag{11}
${\boldsymbol{Q}}$は内力ベクトルになります。
接線剛性マトリクス
有限要素法を解くために、下記の式を満たす節点の変位ベクトル$\boldsymbol{u}$をニュートン法で求めます。
\{\boldsymbol{Q(\boldsymbol{u})}\} = \{\boldsymbol{R}\} \tag{12}
ニュートン法では、接線剛性マトリクス$\boldsymbol{K}=\frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}}$を使用します。
接線剛性マトリクスは 内力ベクトルを変位で線形化することにより求めることができます。
そして、その導出には速度場を導入した線形化手法を用います。
$\boldsymbol{Q}$を速度微分すると、
\frac{d \boldsymbol{Q}}{d t} = \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \frac{d \boldsymbol{u}}{d t} = \boldsymbol{K} \dot{\boldsymbol{u}}\tag{13}
式(13)に$\delta \boldsymbol{u}^T$をかけると、
\delta \boldsymbol{u}^T \frac{d \boldsymbol{Q}}{d t} = \delta \boldsymbol{u}^T \boldsymbol{K} \dot{\boldsymbol{u}} = \delta \boldsymbol{R} \tag{14}
となります。つまり、式(2)を速度微分した式が$\delta \boldsymbol{u}^T \boldsymbol{K} \dot{\boldsymbol{u}}$となります。
式(2)を速度微分すると下記のようになります。
\int_V \dot{(S_{ij}\delta E_{ij} )}dV = \int_V \dot{S_{ij}}\delta E_{ij}dV + \int_V S_{ij} \delta \dot{E}_{ij}dV \tag{15}
ここで、トータルラグランジュ法を簡単に実装する材料構成則としてサンブナン・キルヒホッフモデルがあります。サンブナン・キルヒホッフモデルは下記の式で第2ピオラキルヒホッフ応力を求めることができます。
\begin{align}
&\boldsymbol{S} = \boldsymbol{D} : \boldsymbol{E} \tag{16} \\
&\dot{\boldsymbol{S}} = \boldsymbol{D} : \dot{\boldsymbol{E}} \tag{17}
\end{align}
ここで、$\boldsymbol{D}$は微小変形解析の時と同様のフックの法則のテンソルです。
これを式(15)の右辺第1項に代入すると、
\int_V \dot{S_{ij}}\delta E_{ij}dV = \int_V \boldsymbol{D}_{ijkl} \dot{E_{kl}}\delta E_{ij}dV = \int_V \{\delta \boldsymbol{E}\}^T \bar{\boldsymbol{D}} \{\dot{\boldsymbol{E}}\} dV \tag{18}
となります($\bar{\boldsymbol{D}}$は$\boldsymbol{D}$のvoigt表記)。ここで、$\dot{E_{ij}}$は式(3)を参考に以下のように表すことができるので、
\begin{align}
\dot{E}_{ij} = \frac{1}{2}\{\frac{\partial \dot{u}_i}{\partial X_j} + \frac{\partial \dot{u}_j}{\partial X_i} + (\frac{\partial \dot{u}_k}{\partial X_i} \frac{\partial u_k}{\partial X_j})\ + (\frac{\partial u_k}{\partial X_i} \frac{\partial \dot{u}_k}{\partial X_j})\} \tag{19}
\end{align}
式(8)と同様に
\{ \dot{\boldsymbol{E}}\} = \boldsymbol{B}_{L} \{ \dot{\boldsymbol{u}}\} \tag{20}
式(8)、式(20)を式(18)に代入すると、
\int_V \{\delta \boldsymbol{E}\}^T \bar{\boldsymbol{D}} \{\dot{\boldsymbol{E}}\} dV = \{ \delta \boldsymbol{u}\}^T \int_V \boldsymbol{B}_{L}^T \bar{\boldsymbol{D}} \boldsymbol{B}_{L}dV \{ \dot{\boldsymbol{u}}\} \tag{21}
となります。ここで、$\boldsymbol{K_L}$を以下のように定義します。
\boldsymbol{K_L} = \int_V \boldsymbol{B}_{L}^T \bar{\boldsymbol{D}} \boldsymbol{B}_{L}dV \tag{22}
次に、式(15)の右辺第2項を求めます。$\delta \dot{E}_{ij}$は式(3)を時間微分すれば良いので、
\delta \dot{E}_{ij} = \frac{1}{2}(\frac{\partial \delta u_k}{\partial X_i} \frac{\partial \dot{u}_k}{\partial X_j} + \frac{\partial \dot{u}_k}{\partial X_i} \frac{\partial \delta u_k}{\partial X_j}) \tag{23}
となり、
\int_V S_{ij} \delta \dot{E}_{ij}dV = \frac{1}{2} \int_V (\frac{\partial \delta u_k}{\partial X_i} S_{ij} \frac{\partial \dot{u}_k}{\partial X_j} + \frac{\partial \delta u_k}{\partial X_j} S_{ij} \frac{\partial \dot{u}_k}{\partial X_i}) \tag{24}
と表現することができます。ここで、以下のようにおくと、
\begin{align}
&\boldsymbol{B}_{NL} =
\begin{bmatrix}
\frac{\partial N_1}{\partial X_1} && 0 && 0 && \dots && \frac{\partial N_8}{\partial X_1} && 0 && 0 \\
0 && \frac{\partial N_1}{\partial X_1} && 0 && \dots && 0 && \frac{\partial N_8}{\partial X_1} && 0 \\
0 && 0 && \frac{\partial N_1}{\partial X_1} && \dots && 0 && 0 && \frac{\partial N_8}{\partial X_1} \\
\frac{\partial N_1}{\partial X_2} && 0 && 0 && \dots && \frac{\partial N_8}{\partial X_2} && 0 && 0 \\
0 && \frac{\partial N_1}{\partial X_2} && 0 && \dots && 0 && \frac{\partial N_8}{\partial X_2} && 0 \\
0 && 0 && \frac{\partial N_1}{\partial X_2} && \dots && 0 && 0 && \frac{\partial N_8}{\partial X_2} \\
\frac{\partial N_1}{\partial X_3} && 0 && 0 && \dots && \frac{\partial N_8}{\partial X_3} && 0 && 0 \\
0 && \frac{\partial N_1}{\partial X_3} && 0 && \dots && 0 && \frac{\partial N_8}{\partial X_3} && 0 \\
0 && 0 && \frac{\partial N_1}{\partial X_3} && \dots && 0 && 0 && \frac{\partial N_8}{\partial X_3}
\end{bmatrix}
\\
& \frac{\partial \delta \boldsymbol{u}}{\partial \boldsymbol{X}} = \boldsymbol{B}_{NL} \delta \boldsymbol{u} \\
&\boldsymbol{\bar{S}} =
\begin{bmatrix}
S_{11}\boldsymbol{I} && S_{12}\boldsymbol{I} && S_{13}\boldsymbol{I} \\
S_{21}\boldsymbol{I} && S_{22}\boldsymbol{I} && S_{23}\boldsymbol{I} \\
S_{31}\boldsymbol{I} && S_{32}\boldsymbol{I} && S_{33}\boldsymbol{I}
\end{bmatrix}
\end{align}
式(24)は以下のように表すことができます。
\int_V S_{ij} \delta \dot{E}_{ij}dV = \{ \delta \boldsymbol{u}\}^T \int_V \boldsymbol{B}_{NL}^T \bar{\boldsymbol{S}} \boldsymbol{B}_{NL}dV \{ \dot{\boldsymbol{u}}\} \tag{25}
ここで、$\boldsymbol{K_{NL}}$を以下のように定義します。
\boldsymbol{K_{NL}} = \int_V \boldsymbol{B}_{NL}^T \bar{\boldsymbol{S}} \boldsymbol{B}_{NL} dV \tag{26}
したがって、接線剛性マトリクス$\boldsymbol{K}$は以下のように計算することができます。
\boldsymbol{K} = \boldsymbol{K_L} + \boldsymbol{K_{NL}} \tag{27}
よって、式(11)と式(27)を用いてニュートン法による収束演算をすることができます。
ニュートン法の計算については前回の記事を参照ください。
Pythonで実装する
前回の記事で作成したプログラムをもとに実装します。実装したプログラムはこちらに保存しています。全部のコードを説明すると長くなるので、重要な個所だけ解説します。
サンブナン・キルヒホッフモデルの実装
サンブナン・キルヒホッフモデルの構成則を計算するクラスを実装します。変形勾配テンソル$\boldsymbol{F}$から$\boldsymbol{D}$マトリクスと第2ピオラ・キルヒホッフ応力$\boldsymbol{S}$を計算します。
StVKMaterial.py
import numpy as np
import numpy.linalg as LA
# サンブナン・キルヒホッフモデルの構成則を計算するためのクラス
class StVKMaterial:
# コンストラクタ
# young : パラメータ1
# poisson : パラメータ2
def __init__(self, young, poisson):
# インスタンス変数を定義する
self.young = young
self.poisson = poisson
# Dマトリクスを作成する
# matF : 変形勾配テンソル
def makeDmatrix(self, matF):
tmp = self.young / ((1.0 + self.poisson) * (1.0 - 2.0 * self.poisson))
matD = np.array([[1.0 - self.poisson, self.poisson, self.poisson, 0.0, 0.0, 0.0],
[self.poisson, 1.0 - self.poisson, self.poisson, 0.0, 0.0, 0.0],
[self.poisson, self.poisson, 1.0 - self.poisson, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson), 0.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson), 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poisson)]])
matD = tmp * matD
return matD
# 第2ピオラ・キルヒホッフ応力を作成する
# matF : 変形勾配テンソル
def make2PKStressVector(self, matF):
# グリーン・ラグランジュひずみを計算する
matGLEStrain = 0.5 * (matF.T @ matF - np.eye(3))
vecGLEStrain = np.array([matGLEStrain[0,0], matGLEStrain[1,1], matGLEStrain[2,2],
2.0 * matGLEStrain[0,1], 2.0 * matGLEStrain[1,2], 2.0 * matGLEStrain[0,2]])
# 第2ピオラ・キルヒホッフ応力を計算する
matD = self.makeDmatrix(matF)
vec2PKStress = matD @ vecGLEStrain
return vec2PKStress
要素クラス
makeKetmatrix()で要素剛性マトリクス、makeQeVector()で要素内力ベクトルを計算します。節点の座標が変更されたときはupdate()で要素内の変数を更新します。
C3D8.py
import numpy as np
import numpy.linalg as LA
from ElementOutputData import ElementOutputData
# 6面体8節点要素のクラス
class C3D8:
# コンストラクタ
# no : 要素番号
# nodes : 節点の集合(Node型のリスト)
# matSec : # 材料情報(MaterialSet型)
def __init__(self, no, nodes, matSet):
# インスタンス変数を定義する
self.nodeNum = 8 # 節点の数
self.nodeDof = 3 # 節点の自由度
self.no = no # 要素番号
self.nodes = nodes # 節点の集合(Node型のリスト)
self.density = matSet.density # 密度
self.material = matSet.material # 材料の構成則
self.ipNum = 8 # 積分点の数
self.w1 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # 積分点の重み係数1
self.w2 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # 積分点の重み係数2
self.w3 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] # 積分点の重み係数3
self.ai = np.array([-np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), # 積分点の座標(a,b,c座標系, np.array型のリスト)
-np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0)])
self.bi = np.array([-np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), # 積分点の座標(a,b,c座標系, np.array型のリスト)
-np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0)])
self.ci = np.array([-np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0), # 積分点の座標(a,b,c座標系, np.array型のリスト)
np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0)])
self.matxiyizi = np.array([[self.nodes[0].x, self.nodes[0].y, self.nodes[0].z],
[self.nodes[1].x, self.nodes[1].y, self.nodes[1].z],
[self.nodes[2].x, self.nodes[2].y, self.nodes[2].z],
[self.nodes[3].x, self.nodes[3].y, self.nodes[3].z],
[self.nodes[4].x, self.nodes[4].y, self.nodes[4].z],
[self.nodes[5].x, self.nodes[5].y, self.nodes[5].z],
[self.nodes[6].x, self.nodes[6].y, self.nodes[6].z],
[self.nodes[7].x, self.nodes[7].y, self.nodes[7].z]])
self.incNo = 0 # インクリメントのNo
# 要素内の各行列、ベクトルを初期化する
self.vecDisp = np.zeros(self.nodeNum * self.nodeDof) # 要素内の変位
self.vecStressList = [np.zeros(6) for i in range(self.ipNum)] # 要素内のコーシー応力(np.array型のリスト)
self.misesList = [0.0 for i in range(self.ipNum)] # 要素内のミーゼス応力(リスト)
self.vecGLEStrainList = [np.zeros(6) for i in range(self.ipNum)] # 要素内のグリーン・ラグランジュひずみ(np.array型のリスト)
self.vec2PKStressList = [np.zeros(6) for i in range(self.ipNum)] # 要素内の第2ピオラ・キルヒホッフ応力(np.array型のリスト)
self.centMises = 0.0 # 要素中心のミーゼス応力
self.matdNdabcList = [np.zeros((3,8)) for i in range(self.ipNum)] # dNdabcマトリクス(np.array型のリスト)
self.matdNdxyzList = [np.zeros((3,8)) for i in range(self.ipNum)] # dNdxyzマトリクス(np.array型のリスト)
self.matJList = [np.zeros((3,3)) for i in range(self.ipNum)] # ヤコビ行列(np.array型のリスト)
self.matFList = [np.zeros((3,3)) for i in range(self.ipNum)] # 変形勾配テンソル(np.array型のリスト)
self.matDList = [np.zeros((6,6)) for i in range(self.ipNum)] # Dマトリクス(np.array型のリスト)
self.matBl0List = [np.zeros((6,24)) for i in range(self.ipNum)] # BL0マトリクス(np.array型のリスト)
self.matBl1List = [np.zeros((6,24)) for i in range(self.ipNum)] # BL1マトリクス(np.array型のリスト)
self.matBnlList = [np.zeros((9,24)) for i in range(self.ipNum)] # BNLマトリクス(np.array型のリスト)
self.matSList = [np.zeros((9,9)) for i in range(self.ipNum)] # Sマトリクス(np.array型のリスト)
self.update(np.zeros(self.nodeNum * self.nodeDof), 0)
# 要素接線剛性マトリクスKetを作成する
def makeKetmatrix(self):
# Ketマトリクスをガウス積分で計算する
matKetl = np.zeros([self.nodeDof * self.nodeNum, self.nodeDof * self.nodeNum])
matKetnl = np.zeros([self.nodeDof * self.nodeNum, self.nodeDof * self.nodeNum])
for i in range(self.ipNum):
detMatJ = LA.det(self.matJList[i])
matKetl += self.w1[i] * self.w2[i] * self.w3[i] * \
(self.matBl0List[i] + self.matBl1List[i]).T @ self.matDList[i] @ \
(self.matBl0List[i] + self.matBl1List[i]) * detMatJ
matKetnl += self.w1[i] * self.w2[i] * self.w3[i] * \
self.matBnlList[i].T @ self.matSList[i] @ self.matBnlList[i] * detMatJ
matKet = matKetl + matKetnl
return matKet
# 要素内力ベクトルQeを作成する
def makeQeVector(self):
# 要素内力ベクトルQeを計算する
vecQe = np.zeros(self.nodeDof * self.nodeNum)
for i in range(self.ipNum):
vecQe += self.w1[i] * self.w2[i] * self.w3[i] * \
(self.matBl0List[i] + self.matBl1List[i]).T @ self.vec2PKStressList[i] * LA.det(self.matJList[i])
return vecQe
# 要素内の変数を更新する
# vecDisp : 要素節点の変位ベクトル(np.array型)
# incNo : インクリメントNo
def update(self, vecDisp, incNo):
self.vecDisp = vecDisp
self.incNo = incNo
matuvw = self.vecDisp.reshape(8,3)
for i in range(self.ipNum):
# 各行列、ベクトルを計算する
if incNo == 0:
self.matdNdabcList[i] = self.makedNdabc(self.ai[i], self.bi[i], self.ci[i])
self.matJList[i] = self.makeJmatrix(self.matdNdabcList[i])
self.matdNdxyzList[i] = LA.solve(self.matJList[i], self.matdNdabcList[i])
self.matBl0List[i] = self.makeBl0matrix(self.matdNdxyzList[i])
self.matBnlList[i] = self.makeBnlmatrix(self.matdNdxyzList[i])
matduvwdxyz = self.matdNdxyzList[i] @ matuvw
self.matFList[i] = self.makeFmatrix(matduvwdxyz)
self.vecGLEStrainList[i] = self.makeGLEStrainVector(self.matFList[i])
self.vec2PKStressList[i] = self.make2PKStressVector(self.matFList[i])
self.vecStressList[i] = self.makeStressVector(self.matFList[i], self.vec2PKStressList[i])
self.misesList[i] = self.misesStress(self.vecStressList[i])
self.matDList[i] = self.material.makeDmatrix(self.matFList[i])
self.matBl1List[i] = self.makeBl1matrix(self.matdNdxyzList[i], matduvwdxyz)
self.matSList[i] = self.makeSmatrix(self.vec2PKStressList[i])
# 要素中心値を計算する
if incNo == 0:
matdNdabc = self.makedNdabc(0.0, 0.0, 0.0)
matJ = self.makeJmatrix(matdNdabc)
self.matdNdxyzCent = LA.solve(matJ, matdNdabc)
matduvwdxyz = self.matdNdxyzCent @ matuvw
matF = self.makeFmatrix(matduvwdxyz)
vecCent2PKStress = self.make2PKStressVector(matF)
vecCentStress = self.makeStressVector(matF, vecCent2PKStress)
self.centMises = self.misesStress(vecCentStress)
# 要素の出力データを作成する
def makeOutputData(self):
output = ElementOutputData(self)
return output
# ヤコビ行列を計算する
def makeJmatrix(self, matdNdabc):
# ヤコビ行列を計算する
matJ = matdNdabc @ self.matxiyizi
return matJ
# Bl0マトリクスを作成する
def makeBl0matrix(self, matdNdxyz):
# Bl0マトリクスを計算する
matBl0 = np.zeros((6, self.nodeNum * self.nodeDof))
for i in range(self.nodeNum):
matTmp = np.array([[matdNdxyz[0, i], 0.0, 0.0],
[0.0, matdNdxyz[1, i], 0.0],
[0.0, 0.0, matdNdxyz[2, i]],
[matdNdxyz[1, i], matdNdxyz[0, i], 0.0],
[0.0, matdNdxyz[2, i], matdNdxyz[1, i]],
[matdNdxyz[2, i], 0.0, matdNdxyz[0, i]]])
c = 3 * i
matBl0[:, c:c+3] = matTmp
return matBl0
# Bl1マトリクスを作成する
def makeBl1matrix(self, matdNdxyz, matduvwdxyz):
# Bl1マトリクスを計算する
matBl1 = np.zeros((6, self.nodeNum * self.nodeDof))
for i in range(self.nodeNum):
e11 = matduvwdxyz[0, 0] * matdNdxyz[0, i]
e12 = matduvwdxyz[0, 1] * matdNdxyz[0, i]
e13 = matduvwdxyz[0, 2] * matdNdxyz[0, i]
e21 = matduvwdxyz[1, 0] * matdNdxyz[1, i]
e22 = matduvwdxyz[1, 1] * matdNdxyz[1, i]
e23 = matduvwdxyz[1, 2] * matdNdxyz[1, i]
e31 = matduvwdxyz[2, 0] * matdNdxyz[2, i]
e32 = matduvwdxyz[2, 1] * matdNdxyz[2, i]
e33 = matduvwdxyz[2, 2] * matdNdxyz[2, i]
e41 = matduvwdxyz[1, 0] * matdNdxyz[0, i] + matduvwdxyz[0, 0] * matdNdxyz[1, i]
e42 = matduvwdxyz[1, 1] * matdNdxyz[0, i] + matduvwdxyz[0, 1] * matdNdxyz[1, i]
e43 = matduvwdxyz[1, 2] * matdNdxyz[0, i] + matduvwdxyz[0, 2] * matdNdxyz[1, i]
e51 = matduvwdxyz[1, 0] * matdNdxyz[2, i] + matduvwdxyz[2, 0] * matdNdxyz[1, i]
e52 = matduvwdxyz[1, 1] * matdNdxyz[2, i] + matduvwdxyz[2, 1] * matdNdxyz[1, i]
e53 = matduvwdxyz[1, 2] * matdNdxyz[2, i] + matduvwdxyz[2, 2] * matdNdxyz[1, i]
e61 = matduvwdxyz[2, 0] * matdNdxyz[0, i] + matduvwdxyz[0, 0] * matdNdxyz[2, i]
e62 = matduvwdxyz[2, 1] * matdNdxyz[0, i] + matduvwdxyz[0, 1] * matdNdxyz[2, i]
e63 = matduvwdxyz[2, 2] * matdNdxyz[0, i] + matduvwdxyz[0, 2] * matdNdxyz[2, i]
matTmp = np.array([[e11, e12, e13],
[e21, e22, e23],
[e31, e32, e33],
[e41, e42, e43],
[e51, e52, e53],
[e61, e62, e63]])
c = 3 * i
matBl1[:, c:c+3] = matTmp
return matBl1
# Bnlマトリクスを作成する
def makeBnlmatrix(self, matdNdxyz):
# Bnlマトリクスを計算する
matBnl = np.zeros((9, self.nodeNum * self.nodeDof))
for i in range(self.nodeNum):
matTmp = np.array([[matdNdxyz[0, i], 0.0, 0.0],
[0.0, matdNdxyz[0, i], 0.0],
[0.0, 0.0, matdNdxyz[0, i]],
[matdNdxyz[1, i], 0.0, 0.0],
[0.0, matdNdxyz[1, i], 0.0],
[0.0, 0.0, matdNdxyz[1, i]],
[matdNdxyz[2, i], 0.0, 0.0],
[0.0, matdNdxyz[2, i], 0.0],
[0.0, 0.0, matdNdxyz[2, i]]])
c = 3 * i
matBnl[:, c:c+3] = matTmp
return matBnl
# Sマトリクスを作成する
def makeSmatrix(self, vec2PKStress):
# Sマトリクスを計算する
matS = np.zeros((9,9))
matS11 = np.identity(3) * vec2PKStress[0]
matS22 = np.identity(3) * vec2PKStress[1]
matS33 = np.identity(3) * vec2PKStress[2]
matS12 = np.identity(3) * vec2PKStress[3]
matS23 = np.identity(3) * vec2PKStress[4]
matS31 = np.identity(3) * vec2PKStress[5]
matS[0:3, 0:3] = matS11
matS[0:3, 3:6] = matS12
matS[0:3, 6:9] = matS31
matS[3:6, 0:3] = matS12
matS[3:6, 3:6] = matS22
matS[3:6, 6:9] = matS23
matS[6:9, 0:3] = matS31
matS[6:9, 3:6] = matS23
matS[6:9, 6:9] = matS33
return matS
# dNdabcの行列を計算する
def makedNdabc(self, a, b, c):
# dNi/da, dNi/db, dNi/dcを計算する
dN1da = -0.125 * (1.0 - b) * (1.0 - c)
dN2da = 0.125 * (1.0 - b) * (1.0 - c)
dN3da = 0.125 * (1.0 + b) * (1.0 - c)
dN4da = -0.125 * (1.0 + b) * (1.0 - c)
dN5da = -0.125 * (1.0 - b) * (1.0 + c)
dN6da = 0.125 * (1.0 - b) * (1.0 + c)
dN7da = 0.125 * (1.0 + b) * (1.0 + c)
dN8da = -0.125 * (1.0 + b) * (1.0 + c)
dN1db = -0.125 * (1.0 - a) * (1.0 - c)
dN2db = -0.125 * (1.0 + a) * (1.0 - c)
dN3db = 0.125 * (1.0 + a) * (1.0 - c)
dN4db = 0.125 * (1.0 - a) * (1.0 - c)
dN5db = -0.125 * (1.0 - a) * (1.0 + c)
dN6db = -0.125 * (1.0 + a) * (1.0 + c)
dN7db = 0.125 * (1.0 + a) * (1.0 + c)
dN8db = 0.125 * (1.0 - a) * (1.0 + c)
dN1dc = -0.125 * (1.0 - a) * (1.0 - b)
dN2dc = -0.125 * (1.0 + a) * (1.0 - b)
dN3dc = -0.125 * (1.0 + a) * (1.0 + b)
dN4dc = -0.125 * (1.0 - a) * (1.0 + b)
dN5dc = 0.125 * (1.0 - a) * (1.0 - b)
dN6dc = 0.125 * (1.0 + a) * (1.0 - b)
dN7dc = 0.125 * (1.0 + a) * (1.0 + b)
dN8dc = 0.125 * (1.0 - a) * (1.0 + b)
# dNdabcを計算する
dNdabc = np.array([[dN1da, dN2da, dN3da, dN4da, dN5da, dN6da, dN7da, dN8da],
[dN1db, dN2db, dN3db, dN4db, dN5db, dN6db, dN7db, dN8db],
[dN1dc, dN2dc, dN3dc, dN4dc, dN5dc, dN6dc, dN7dc, dN8dc]])
return dNdabc
# 変形勾配テンソルを作成する
def makeFmatrix(self, matduvwdxyz):
# 変形勾配テンソルを計算する
matF = np.eye(3) + matduvwdxyz.T
return matF
# グリーン・ラグランジュひずみベクトルを作成する
def makeGLEStrainVector(self, matF):
# グリーン・ラグランジュひずみを計算する
matGLEStrain = 0.5 * (matF.T @ matF - np.eye(3))
vecGLEStrain = np.array([matGLEStrain[0,0], matGLEStrain[1,1], matGLEStrain[2,2], 2.0 * matGLEStrain[0,1], 2.0 * matGLEStrain[1,2], 2.0 * matGLEStrain[0,2]])
return vecGLEStrain
# 第2ピオラ・キルヒホッフ応力ベクトルを作成する
def make2PKStressVector(self, matF):
# グリーン・ラグランジュひずみを計算する
vec2PKStress = self.material.make2PKStressVector(matF)
return vec2PKStress
# コーシー応力ベクトルを作成する
def makeStressVector(self, matF, vec2PKStress):
# コーシー応力ベクトルを計算する
S11 = vec2PKStress[0]
S22 = vec2PKStress[1]
S33 = vec2PKStress[2]
S12 = vec2PKStress[3]
S23 = vec2PKStress[4]
S13 = vec2PKStress[5]
mat2PKStress = np.array([[S11, S12, S13],
[S12, S22, S23],
[S13, S23, S33]])
matStress = matF @ mat2PKStress @ matF.T / LA.det(matF)
vecStress = np.array([matStress[0, 0], matStress[1, 1], matStress[2, 2], matStress[0, 1], matStress[1, 2], matStress[0,2]])
return vecStress
# ミーゼス応力を計算する
# vecStress : 応力ベクトル(np.array型)
def misesStress(self, vecStress):
tmp1 = np.square(vecStress[0] - vecStress[1]) + np.square(vecStress[1] - vecStress[2]) + np.square(vecStress[2] - vecStress[0])
tmp2 = 6.0 * (np.square(vecStress[3]) + np.square(vecStress[4]) + np.square(vecStress[5]))
mises = np.sqrt(0.5 * (tmp1 + tmp2))
return mises
# 体積を求める
def getVolume(self):
# ガウス積分で体積を計算する
volume = 0
for i in range(self.ipNum):
volume += self.w1[i] * self.w2[i] * self.w3[i] * LA.det(self.matJList[i])
return volume
FEM解析を行うクラス
FEM解析を行うクラスです。前回の記事のコードと基本的には一緒ですが、scipyによる高速化、(ParaViewで確認用の)VTKファイル出力機能を追加しています。
import numpy as np
import numpy.linalg as LA
from scipy.sparse import lil_matrix, coo_matrix
from scipy.sparse.linalg import spsolve
from OutputVTK import OutputVTK
class FEM:
# コンストラクタ
# nodes : 節点は1から始まる順番で並んでいる前提(Node型のリスト)
# elements : 要素は種類ごとにソートされている前提(C3D8型のリスト)
# bound : 境界条件(d2Boundary型)
# incNum : インクリメント数
def __init__(self, nodes, elements, bound, incNum):
# インスタンス変数を定義する
self.nodeDof = 3 # 節点の自由度
self.nodes = nodes # 節点は1から始まる順番で並んでいる前提(リスト)
self.elements = elements # 要素は種類ごとにソートされている前提(リスト)
self.bound = bound # 境界条件(Boundary型)
self.incNum = incNum # インクリメント数
self.itrNum = 100 # ニュートン法のイテレータの上限
self.rn = 0.005 # ニュートン法の残差力の収束判定のパラメータ
#self.cn = 0.01 # ニュートン法の変位の収束判定のパラメータ
self.cn = 0.001 # ニュートン法の変位の収束判定のパラメータ
# 陰解法で解析を行う
def impAnalysis(self):
self.vecDispList = [] # インクリメント毎の変位ベクトルのリスト(np.array型のリスト)
self.vecRFList = [] # インクリメント毎の反力ベクトルのリスト(np.array型のリスト)
self.elemOutputDataList = [] # インクリメント毎の要素出力のリスト(makeOutputData型のリストのリスト)
# 荷重をインクリメント毎に分割する
vecfList = []
for i in range(self.incNum):
vecfList.append(self.makeForceVector() * (i + 1) / self.incNum)
# ステータス表示用
print("------------")
print("Status")
print("------------")
print("INC".rjust(6) + "ITR".rjust(6))
# ニュートン法により変位を求める
vecDisp = np.zeros(len(self.nodes) * self.nodeDof) # 全節点の変位ベクトル
vecR = np.zeros(len(self.nodes) * self.nodeDof) # 残差力ベクトル
for i in range(self.incNum):
vecDispFirst = vecDisp.copy() # 初期の全節点の変位ベクトル
vecf = vecfList[i] # i+1番インクリメントの荷重
vecBoundDisp = self.bound.makeDispVector()
# 境界条件を考慮しないインクリメント初期の残差力ベクトルRを作成する
if i == 0:
vecR = vecfList[0]
else:
vecR = vecfList[i] - vecfList[i - 1]
# 収束演算を行う
matKt = self.makeKtmatrix() # 接線剛性マトリクスKt(CSR)を作成する
matKtc, vecRc = self.setBoundCondition(matKt, vecR, vecBoundDisp, vecDisp) # 境界条件を考慮したKtcマトリクス、Rcベクトルを作成する
for j in range(self.itrNum):
# 変位ベクトルを計算する
vecd = spsolve(matKtc, vecRc)
if not np.isfinite(vecd).all():
raise ValueError("有限要素法の計算に失敗しました。")
vecDisp += vecd
# 全ての要素内の変数を更新する
self.updateElements(vecDisp, i+1)
# 新たな残差力ベクトルRを求める
vecQ = self.makeVecQ()
vecR = vecf - vecQ
# 新たな境界条件を考慮したRcベクトルを作成する
matKt = self.makeKtmatrix()
matKtc, vecRc = self.setBoundCondition(matKt, vecR, vecBoundDisp, vecDisp)
# 時間平均力を計算する
aveForce = 0.0
cnt = len(self.nodes) * self.nodeDof
for k in range(len(vecQ)):
aveForce += np.abs(vecQ[k])
for k in range(len(vecf)):
if not vecf[k] == 0.0:
aveForce += np.abs(vecf[k])
cnt+= 1
aveForce = aveForce / cnt
# 収束判定を行う
if np.allclose(vecRc, 0.0):
break
if np.isclose(LA.norm(vecd), 0.0):
break
dispRate = LA.norm(vecd) / LA.norm(vecDisp - vecDispFirst)
ResiForceRate = np.abs(vecRc).max() / aveForce
if dispRate < self.cn and ResiForceRate < self.rn:
break
# 収束できなかった場合
if j == (self.itrNum - 1):
raise ValueError("収束しませんでした。")
# インクリメントの最終的な変位べクトルを格納する
self.vecDispList.append(vecDisp.copy())
# 変位ベクトルから要素の出力データを計算する
elemOutputDatas = []
for elem in self.elements:
elemOutputData = elem.makeOutputData()
elemOutputDatas.append(elemOutputData)
self.elemOutputDataList.append(elemOutputDatas)
# 節点反力を計算する
vecRF = np.array(vecQ - vecf).flatten()
# インクリメントの最終的な節点反力を格納する
self.vecRFList.append(vecRF)
# ステータスを表示する
print(str(i+1).rjust(6) + str(j+1).rjust(6))
# 完了したことを表示する
print("The analysis has completed")
# 接線剛性マトリクスKtを作成する(CSR形式)
def makeKtmatrix(self):
nDof = len(self.nodes) * self.nodeDof
elemNum = len(self.elements)
nnz = elemNum * 24 * 24
rows = np.empty(nnz, dtype=np.int32)
cols = np.empty(nnz, dtype=np.int32)
data = np.empty(nnz)
p = 0
for elem in self.elements:
# ketマトリクスを計算する
matKet = elem.makeKetmatrix()
dofs = np.empty(24, dtype=np.int32)
for i in range(8):
base = (elem.nodes[i].no - 1) * 3
dofs[3 * i] = base
dofs[3 * i + 1] = base + 1
dofs[3 * i + 2] = base + 2
for i in range(24):
for j in range(24):
rows[p] = dofs[i]
cols[p] = dofs[j]
data[p] = matKet[i,j]
p += 1
return coo_matrix((data,(rows,cols)),shape=(nDof,nDof)).tocsr()
# 節点に負荷する荷重、等価節点力を考慮した荷重ベクトルを作成する
def makeForceVector(self):
# 節点に負荷する荷重ベクトルを作成する
vecCondiForce = self.bound.makeForceVector()
vecf = vecCondiForce
return vecf
# Kマトリクス、荷重ベクトルに境界条件を考慮する
# matKt : 接線剛性マトリクス
# vecR : 残差力ベクトル
# vecBoundDisp : 節点の境界条件の変位ベクトル
# vecDisp : 全節点の変位ベクトル(np.array型)
def setBoundCondition(self, matKt, vecR, vecBoundDisp, vecDisp):
matKtc = matKt.tolil().copy()
vecRc = np.copy(vecR)
# 単点拘束条件を考慮したKマトリクス、荷重ベクトルを作成する
for i in range(len(vecBoundDisp)):
if not vecBoundDisp[i] == None:
# Kマトリクスからi列を抽出する
vecx = np.array(matKt[:, i].todense()).flatten()
# 変位ベクトルi列の影響を荷重ベクトルに適用する
vecRc = vecRc - (vecBoundDisp[i] - vecDisp[i]) * vecx
# Kマトリクスのi行、i列を全て0にし、i行i列の値を1にする
matKtc[:, i] = 0.0
matKtc[i, :] = 0.0
matKtc[i, i] = 1.0
matKtc = matKtc.tocsr()
for i in range(len(vecBoundDisp)):
if not vecBoundDisp[i] == None:
vecRc[i] = vecBoundDisp[i] - vecDisp[i]
return matKtc, vecRc
# 全ての要素内の変数を更新する
# vecDisp : 全節点の変位ベクトル(np.array型)
# incNo : インクリメントの番号
def updateElements(self, vecDisp, incNo):
for elem in self.elements:
vecElemDisp = np.zeros(len(elem.nodes) * self.nodeDof)
for i in range(len(elem.nodes)):
for j in range(elem.nodeDof):
vecElemDisp[i * elem.nodeDof + j] = vecDisp[(elem.nodes[i].no - 1) * self.nodeDof + j]
elem.update(vecElemDisp, incNo)
# 内力ベクトルを作成する
def makeVecQ(self):
vecQ = np.zeros(len(self.nodes) * self.nodeDof)
for elem in self.elements:
vecQe = elem.makeQeVector()
for k in range(len(elem.nodes)):
for l in range(elem.nodeDof):
vecQ[(elem.nodes[k].no - 1) * self.nodeDof + l] += vecQe[k * elem.nodeDof + l]
return vecQ
# 解析結果をテキストファイルに出力する
def outputTxt(self, filePath):
# ファイルを作成し、開く
f = open(filePath + ".txt", 'w')
# 出力する文字の情報を定義する
columNum = 20
floatDigits = ".10g"
# 入力データのタイトルを書きこむ
f.write("*********************************\n")
f.write("* Input Data *\n")
f.write("*********************************\n")
f.write("\n")
# 節点情報を出力する
f.write("***** Node Data ******\n")
f.write("No".rjust(columNum) + "X".rjust(columNum) + "Y".rjust(columNum) + "Z".rjust(columNum) + "\n")
f.write("-" * columNum * 4 + "\n")
for node in self.nodes:
strNo = str(node.no).rjust(columNum)
strX = str(format(node.x, floatDigits).rjust(columNum))
strY = str(format(node.y, floatDigits).rjust(columNum))
strZ = str(format(node.z, floatDigits).rjust(columNum))
f.write(strNo + strX + strY + strZ + "\n")
f.write("\n")
# 要素情報を出力する
nodeNoColumNum = 36
f.write("***** Element Data ******\n")
f.write("No".rjust(columNum) + "Type".rjust(columNum) + "Node No".rjust(nodeNoColumNum) + "Density".rjust(columNum) + "\n")
f.write("-" * columNum * 7 + "-" * nodeNoColumNum + "\n")
for elem in self.elements:
strNo = str(elem.no).rjust(columNum)
strType = str(elem.__class__.__name__ ).rjust(columNum)
strNodeNo = ""
for node in elem.nodes:
strNodeNo += " " + str(node.no)
strNodeNo = strNodeNo.rjust(nodeNoColumNum)
strDensity = "None".rjust(columNum)
if not elem.density is None:
strDensity = str(format(elem.density, floatDigits).rjust(columNum))
f.write(strNo + strType + strNodeNo + strDensity + "\n")
f.write("\n")
# 単点拘束情報を出力する
f.write("***** SPC Constraint Data ******\n")
f.write("NodeNo".rjust(columNum) + "X Displacement".rjust(columNum) + "Y Displacement".rjust(columNum) +
"Z Displacement".rjust(columNum) + "\n")
f.write("-" * columNum * 4 + "\n")
vecBoundDisp = self.bound.makeDispVector()
for i in range(self.bound.nodeNum):
strFlg = False
for j in range(self.bound.nodeDof):
if not vecBoundDisp[i * self.bound.nodeDof + j] is None:
strFlg = True
if strFlg == True:
strNo = str(i + 1).rjust(columNum)
strXDisp = "None".rjust(columNum)
if not vecBoundDisp[i * self.bound.nodeDof] is None:
strXDisp = str(format(vecBoundDisp[i * self.bound.nodeDof], floatDigits).rjust(columNum))
strYDisp = "None".rjust(columNum)
if not vecBoundDisp[i * self.bound.nodeDof + 1] is None:
strYDisp = str(format(vecBoundDisp[i * self.bound.nodeDof + 1], floatDigits).rjust(columNum))
strZDisp = "None".rjust(columNum)
if not vecBoundDisp[i * self.bound.nodeDof + 2] is None:
strZDisp = str(format(vecBoundDisp[i * self.bound.nodeDof + 2], floatDigits).rjust(columNum))
f.write(strNo + strXDisp + strYDisp + strZDisp + "\n")
f.write("\n")
# 荷重条件を出力する(等価節点力も含む)
f.write("***** Nodal Force Data ******\n")
f.write("NodeNo".rjust(columNum) + "X Force".rjust(columNum) + "Y Force".rjust(columNum) +
"Z Force".rjust(columNum) + "\n")
f.write("-" * columNum * 4 + "\n")
vecf = self.makeForceVector()
for i in range(len(self.nodes)):
strFlg = False
for j in range(self.bound.nodeDof):
if not vecf[i * self.bound.nodeDof + j] == 0.0:
strFlg = True
if strFlg == True:
strNo = str(i + 1).rjust(columNum)
strXForce = str(format(vecf[i * self.bound.nodeDof], floatDigits).rjust(columNum))
strYForce = str(format(vecf[i * self.bound.nodeDof + 1], floatDigits).rjust(columNum))
strZForce = str(format(vecf[i * self.bound.nodeDof + 2], floatDigits).rjust(columNum))
f.write(strNo + strXForce + strYForce + strZForce + "\n")
f.write("\n")
# 結果データのタイトルを書きこむ
f.write("**********************************\n")
f.write("* Result Data *\n")
f.write("**********************************\n")
f.write("\n")
for i in range(self.incNum):
f.write("*Increment " + str(i + 1) + "\n")
f.write("\n")
# 変位のデータを出力する
f.write("***** Displacement Data ******\n")
f.write("NodeNo".rjust(columNum) + "Magnitude".rjust(columNum) + "X Displacement".rjust(columNum) +
"Y Displacement".rjust(columNum) + "Z Displacement".rjust(columNum) + "\n")
f.write("-" * columNum * 5 + "\n")
for j in range(len(self.nodes)):
strNo = str(j + 1).rjust(columNum)
vecDisp = self.vecDispList[i]
mag = np.linalg.norm(np.array((vecDisp[self.nodeDof * j], vecDisp[self.nodeDof * j + 1], vecDisp[self.nodeDof * j + 2])))
strMag = str(format(mag, floatDigits).rjust(columNum))
strXDisp = str(format(vecDisp[self.nodeDof * j], floatDigits).rjust(columNum))
strYDisp = str(format(vecDisp[self.nodeDof * j + 1], floatDigits).rjust(columNum))
strZDisp = str(format(vecDisp[self.nodeDof * j + 2], floatDigits).rjust(columNum))
f.write(strNo + strMag + strXDisp + strYDisp + strZDisp + "\n")
f.write("\n")
# 応力データを出力する
f.write("***** Stress Data ******\n")
f.write("Element No".rjust(columNum) + "Integral No".rjust(columNum) + "Stress XX".rjust(columNum) + "Stress YY".rjust(columNum) +
"Stress ZZ".rjust(columNum) + "Stress XY".rjust(columNum) + "Stress XZ".rjust(columNum) + "Stress YZ".rjust(columNum) +
"Mises".rjust(columNum) + "\n")
f.write("-" * columNum * 9 + "\n")
for elemOutputData in self.elemOutputDataList[i]:
elem = elemOutputData.element
strElemNo = str(elem.no).rjust(columNum)
for j in range(elem.ipNum):
strIntNo = str(j + 1).rjust(columNum)
strStressXX = str(format(elemOutputData.vecIpStressList[j][0], floatDigits).rjust(columNum))
strStressYY = str(format(elemOutputData.vecIpStressList[j][1], floatDigits).rjust(columNum))
strStressZZ = str(format(elemOutputData.vecIpStressList[j][2], floatDigits).rjust(columNum))
strStressXY = str(format(elemOutputData.vecIpStressList[j][3], floatDigits).rjust(columNum))
strStressXZ = str(format(elemOutputData.vecIpStressList[j][4], floatDigits).rjust(columNum))
strStressYZ = str(format(elemOutputData.vecIpStressList[j][5], floatDigits).rjust(columNum))
strMises = str(format(elemOutputData.ipMiseses[j], floatDigits).rjust(columNum))
f.write(strElemNo + strIntNo + strStressXX + strStressYY + strStressZZ +
strStressXY + strStressXZ + strStressYZ + strMises + "\n")
f.write("\n")
# 第2ピオラ・キルヒホッフ応力データを出力する
f.write("***** 2PK Stress Data ******\n")
f.write("Element No".rjust(columNum) + "Integral No".rjust(columNum) + "2PKStress XX".rjust(columNum) + "2PKStress YY".rjust(columNum) +
"2PKStress ZZ".rjust(columNum) + "2PKStress XY".rjust(columNum) + "2PKStress XZ".rjust(columNum) + "2PKStress YZ".rjust(columNum) + "\n")
f.write("-" * columNum * 9 + "\n")
for elemOutputData in self.elemOutputDataList[i]:
elem = elemOutputData.element
strElemNo = str(elem.no).rjust(columNum)
for j in range(elem.ipNum):
strIntNo = str(j + 1).rjust(columNum)
str2PKStressXX = str(format(elemOutputData.vecIp2PKStressList[j][0], floatDigits).rjust(columNum))
str2PKStressYY = str(format(elemOutputData.vecIp2PKStressList[j][1], floatDigits).rjust(columNum))
str2PKStressZZ = str(format(elemOutputData.vecIp2PKStressList[j][2], floatDigits).rjust(columNum))
str2PKStressXY = str(format(elemOutputData.vecIp2PKStressList[j][3], floatDigits).rjust(columNum))
str2PKStressXZ = str(format(elemOutputData.vecIp2PKStressList[j][4], floatDigits).rjust(columNum))
str2PKStressYZ = str(format(elemOutputData.vecIp2PKStressList[j][5], floatDigits).rjust(columNum))
f.write(strElemNo + strIntNo + str2PKStressXX + str2PKStressYY + str2PKStressZZ +
str2PKStressXY + str2PKStressXZ + str2PKStressYZ + "\n")
f.write("\n")
# グリーン・ラグランジュひずみデータを出力する
f.write("***** Green-Lagrange Strain Data ******\n")
f.write("Element No".rjust(columNum) + "Integral No".rjust(columNum) + "GLStrain XX".rjust(columNum) + "GLStrain YY".rjust(columNum) +
"GLStrain ZZ".rjust(columNum) + "GLStrain XY".rjust(columNum) + "GLStrain XZ".rjust(columNum) + "GLStrain YZ".rjust(columNum) + "\n")
f.write("-" * (columNum * 8 + 30) + "\n")
for elemOutputData in self.elemOutputDataList[i]:
elem = elemOutputData.element
strElemNo = str(elem.no).rjust(columNum)
for j in range(elem.ipNum):
strIntNo = str(j + 1).rjust(columNum)
strGLStrainXX = str(format(elemOutputData.vecIpGLEStrainList[j][0], floatDigits).rjust(columNum))
strGLStrainYY = str(format(elemOutputData.vecIpGLEStrainList[j][1], floatDigits).rjust(columNum))
strGLStrainZZ = str(format(elemOutputData.vecIpGLEStrainList[j][2], floatDigits).rjust(columNum))
strGLStrainXY = str(format(elemOutputData.vecIpGLEStrainList[j][3], floatDigits).rjust(columNum))
strGLStrainXZ = str(format(elemOutputData.vecIpGLEStrainList[j][4], floatDigits).rjust(columNum))
strGLStrainYZ = str(format(elemOutputData.vecIpGLEStrainList[j][5], floatDigits).rjust(columNum))
f.write(strElemNo + strIntNo + strGLStrainXX + strGLStrainYY + strGLStrainZZ +
strGLStrainXY + strGLStrainXZ + strGLStrainYZ + "\n")
f.write("\n")
# 反力のデータを出力する
f.write("***** Reaction Force Data ******\n")
f.write("NodeNo".rjust(columNum) + "Magnitude".rjust(columNum) + "X Force".rjust(columNum) + "Y Force".rjust(columNum) +
"Z Force".rjust(columNum) + "\n")
f.write("-" * columNum * 5 + "\n")
for j in range(len(self.nodes)):
strNo = str(j + 1).rjust(columNum)
vecRF = self.vecRFList[i]
mag = np.linalg.norm(np.array((vecRF[self.nodeDof * j], vecRF[self.nodeDof * j + 1], vecRF[self.nodeDof * j + 2])))
strMag = str(format(mag, floatDigits).rjust(columNum))
strXForce = str(format(vecRF[self.nodeDof * j], floatDigits).rjust(columNum))
strYForce = str(format(vecRF[self.nodeDof * j + 1], floatDigits).rjust(columNum))
strZForce = str(format(vecRF[self.nodeDof * j + 2], floatDigits).rjust(columNum))
f.write(strNo + strMag + strXForce + strYForce + strZForce + "\n")
f.write("\n")
# ファイルを閉じる
f.close()
# 解析結果をvtkファイルに出力する
def outputVTK(self, filePath):
outputVTK = OutputVTK(self.nodes, self.elements, self.vecDispList, self.vecRFList, self.elemOutputDataList, self.incNum)
outputVTK.output(filePath)
サンプルを解く
以下の図の片持ち梁のサンプルを解き、理論値のたわみ量を比較します。
ここで、L=500mm、h=10mm、b=10mmでFは10Nと5000Nとした2ケース試してみます。
今回比較している梁理論のたわみ式は微小変形理論に基づく解であり、今回実装したサンブナン・キルヒホッフモデルの厳密解ではありません。そのため、大変形領域では一致しません。一方で荷重が十分小さい場合は有限変形理論も微小変形理論へ漸近するため、実装検証として利用することができます。
サンプルを解くソースコードを下記に示します。
main.py
import copy
import time
from Node import Node
from MaterialSet import MaterialSet
from StVKMaterial import StVKMaterial
from Boundary import Boundary
from FEM import FEM
from C3D8 import C3D8
# メインの処理
def main():
# ============================
# メッシュ設定
# ============================
nx = 2 # x方向要素数
ny = 2 # y方向要素数
nz = 200 # z方向要素数
lx = 10.0
ly = 10.0
lz = 500.0
dx = lx / nx
dy = ly / ny
dz = lz / nz
# ============================
# 節点生成
# ============================
node_dict = {}
node_id = 1
for ix in range(nx + 1):
x = -lx / 2 + ix * dx
for iy in range(ny + 1):
y = -ly / 2 + iy * dy
for iz in range(nz + 1):
z = lz - iz * dz
node_dict[node_id] = Node(node_id,x,y,z)
node_id += 1
nodes = [node_dict[i] for i in range(1, node_id)]
# ============================
# 節点番号取得関数
# ============================
nodes_per_y = ny + 1
nodes_per_z = nz + 1
def nid(ix, iy, iz):
return (ix * nodes_per_y * nodes_per_z + iy * nodes_per_z + iz+ 1)
# ============================
# 材料
# ============================
young = 210000.0
poisson = 0.0
mat = StVKMaterial(young,poisson)
matSet = MaterialSet("None",material=mat)
# ============================
# 要素生成
# ============================
elems = []
elem_id = 1
for ix in range(nx):
for iy in range(ny):
for iz in range(nz):
n1 = node_dict[nid(ix , iy , iz+1)]
n2 = node_dict[nid(ix+1, iy , iz+1)]
n3 = node_dict[nid(ix+1, iy+1, iz+1)]
n4 = node_dict[nid(ix , iy+1, iz+1)]
n5 = node_dict[nid(ix , iy , iz)]
n6 = node_dict[nid(ix+1, iy , iz)]
n7 = node_dict[nid(ix+1, iy+1, iz)]
n8 = node_dict[nid(ix , iy+1, iz)]
elems.append(C3D8(elem_id,[n1, n2, n3, n4,n5, n6, n7, n8],matSet))
elem_id += 1
# ============================
# 境界条件
# ============================
bound = Boundary(len(nodes))
# 固定端(z=0)
for ix in range(nx + 1):
for iy in range(ny + 1):
node = node_dict[nid(ix, iy, nz)]
bound.addSPC(node.no,0.0,0.0,0.0)
# 自由端中央
bound1 = copy.deepcopy(bound)
bound2 = copy.deepcopy(bound)
tip = node_dict[nid(nx // 2,ny // 2,0)]
bound1.addForce(tip.no,10.0,0.0,0.0)
bound2.addForce(tip.no,5000.0,0.0,0.0)
# ============================
# FEM解析
# ============================
fem1 = FEM(nodes,elems,bound1,10)
start = time.perf_counter() #計測開始
fem1.impAnalysis()
end = time.perf_counter() #計測終了
print("解析時間 : " + '{:.2f}'.format((end-start)/60) + "分")
# ============================
# 結果出力
# ============================
fem1.outputTxt("test_10N")
fem1.outputVTK("test_10N")
# ============================
# FEM解析
# ============================
fem2 = FEM(nodes,elems,bound2,10)
start = time.perf_counter() #計測開始
fem2.impAnalysis()
end = time.perf_counter() #計測終了
print("解析時間 : " + '{:.2f}'.format((end-start)/60) + "分")
# ============================
# 結果出力
# ============================
fem2.outputTxt("test_5000N")
fem2.outputVTK("test_5000N")
if __name__ == '__main__':
main()
解析モデル
断面は2×2、長さは200分割でメッシュを作成しました。せん断ロッキングを抑えるために長さ方向に細かめのメッシュ分割にしています。
解析結果
F=10[N]では、以下の結果が得られました。
たわみ量
| 理論解 | FEM | 誤差率 |
|---|---|---|
| 2.381[mm] | 2.309[mm] | 3.02[%] |
材料の構成則が理論解と違いますが、微小な変形なら概ね一致はします。
F=5000[N]とした場合は、大変形しますが、それっぽい結果は得られます。
(ここまで変形すると微小変形解析では変形挙動がおかしくなります。)
まとめ
本記事では
・Green-Lagrangeひずみ
・第2Piola-Kirchhoff応力
・Total Lagrangian法
・Saint-Venant Kirchhoffモデル
を用いて有限変形FEMを実装し、その検証をしてみました。
参考文献
[1]非線形有限要素法特論
[2]非線形有限要素法の基礎と応用