【python】有限要素法(FEM)をゼロから実装する2_関係式の整理

はじめに

前回から有限要素法による解析プログラムを自分で作成してみています。
前回はメッシュを作成するところまでやったので、今回はプログラムを書くにあたって重要になる材料力学上の関係式を整理します。
重要な式だけ抜き出しましたが、数式が多くて今回はコードなしです。

はじめに
関係式の整理
終わりに
参考文献

関係式の整理

まずは静弾性問題における要素内の関係式を整理します。

応力とひずみの関係式

要素内の任意の点に働く応力ベクトルを $\boldsymbol{\sigma}$ 、ひずみベクトルを $\boldsymbol{\epsilon}$ とすると以下の関係になります。

$\boldsymbol{\sigma} = \boldsymbol{D\epsilon}$

ここで $\boldsymbol{D}$ は応力-ひずみマトリクス(Dマトリクス)と呼ばれ、材料特性によって決まる値です。

今回の2次元問題の場合は、2次元近似の方法の違いで以下の2種類があります。

平面応力(plane-stress): Ｚ方向の応力をゼロと仮定したもので主に薄い板の解析に使用される
平面ひずみ(plane-strain): Ｚ方向のひずみをゼロと仮定したものでZ方向に大きいものに使用される(棒とかトンネルとかを輪切りにして解析するイメージ)

今回は薄い板の解析を対象にしたいと思うので平面応力状態とします。その場合、応力とひずみの関係は以下のようになります。

$\begin{bmatrix} \sigma_{x} \\ \sigma_{y} \\ \tau_{xy} \end{bmatrix} = \dfrac{E}{(1-\nu)(1+\nu)}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}\begin{bmatrix} \epsilon_{x} \\ \epsilon_{y} \\ \gamma_{xy} \end{bmatrix}$

ここで $E$ はヤング率、 $\nu$ はポアソン比です。よって $\boldsymbol{D}$ は以下で求められます。

$\boldsymbol{D} = \dfrac{E}{(1-\nu)(1+\nu)}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix}$

ひずみと変位の関係式

要素内の任意の点におけるひずみと変位ベクトル $\boldsymbol{u}$ は以下の関係式で表すことができます。
$\boldsymbol{\epsilon} = \boldsymbol{Lu}$

ここで $\boldsymbol{L}$ は以下と置いています。

$\boldsymbol{L} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} \end{bmatrix}$

変位ベクトル $\boldsymbol{u}$ は要素内の任意の位置の変位です。

仮想仕事の原理の式

要素内で微小な変位 $\boldsymbol{\delta u}$ を仮定した場合、その変異によってなされる仕事は釣り合うという関係式です。ここで要素に働く体積力(重力とか)はないとすると下記のようにあらわすことができます。

$\int \boldsymbol{\delta u}^{T}\boldsymbol{P}dS = \iint \boldsymbol{\delta\epsilon}^{T}\boldsymbol{\sigma}hdxdy$

ここで、 $h$ は要素の厚み、 $\boldsymbol{P}$ は要素に働く表面力 $\boldsymbol{\delta\epsilon}$ は仮想ひずみベクトルです。
先の2式を使って式変形すると、以下の形になります。

$\boldsymbol{\delta u}^{T}\boldsymbol{f_{e}} = h \iint \boldsymbol{\delta u}^{T}\boldsymbol{L}^{T}\boldsymbol{DLu}dxdy$

$\boldsymbol{f_{e}}$ は等価節点力といい、各節点に働く荷重をまとめたベクトルです。
要素内でdxdyの2重積分がありますが、この積分を簡単に計算するために全部の要素が同じ整った形状になるような座標変換(正規化)を行います。

形状関数の導入

具体的には以下のように、 $(x, y)$ の座標系から、各要素が以下の正方形になるような $(\xi, \eta)$ 座標系に変換します。

この座標系 $(\xi, \eta)$ で考えることで積分が非常に簡単になります。

座標系 $(x, y)$ における要素の節点の座標をまとめた節点座標ベクトルを $\boldsymbol{q}$ とします。要素内の任意の点の座標ベクトル $\boldsymbol{p}$ を $(\xi, \eta)$ 座標系で表した場合、以下のようになります。

$\boldsymbol{p}(\xi, \eta) = \begin{bmatrix} x(\xi, \eta) \\ y(\xi, \eta) \\ \end{bmatrix} = \boldsymbol{Nq}(x, y)$

ここで各ノードの座標( $(x, y)$ 座標系)を $(x_{i}, y_{i})$ とすると節点座標ベクトル $\boldsymbol{x}$ 、形状関数 $\boldsymbol{N}$ は以下のようになります(4角形1次要素の場合)。

$\boldsymbol{x} = \begin{bmatrix} x_{1} \\ y_{1} \\ x_{2} \\ y_{2} \\ x_{3} \\ y_{3} \\ x_{4} \\ y_{4} \\ \end{bmatrix}, \boldsymbol{N} = \begin{bmatrix} N_{1} & 0 & N_{2} & 0 & N_{3} & 0 & N_{4} & 0 \\ 0 & N_{1} & 0 & N_{2} & 0 & N_{3} & 0 & N_{4} \\ \end{bmatrix}$

$N_{i}$ は以下になります。

$N_{1} = \dfrac{1}{4}(1 - \xi)(1 - \eta)$
$N_{2} = \dfrac{1}{4}(1 + \xi)(1 - \eta)$
$N_{3} = \dfrac{1}{4}(1 + \xi)(1 + \eta)$
$N_{4} = \dfrac{1}{4}(1 - \xi)(1 + \eta)$

先程4角形1次要素という話をしましたが、4角形1次要素というのは、4つのノードで1つの要素を表す要素です。計算が簡単ですが、曲げとかのモードには精度が低いです。
4角形2次要素もあり、これは4つのノードにそれぞれの中点4つも加えた8つのノードで表します。こちらは精度がいいのでよく使われます。

ただ、今回は簡単なので4角形1次要素で行きます。2次要素に関しては形状関数の形とか色々でまわっているので調べてみてください。

形状関数の変位への適用

先ほどは要素の座標において形状関数を用いて $(\xi, \eta)$ 座標系に変更しましたが、要素内の任意の点の変位も同じ式を適用できます。
つまり、各節点の変位を用いて要素内任意の点の $(\xi, \eta)$ 座標系における変位ベクトルを $\boldsymbol{u}$ とした場合に、下記の式で表すことができます。

$\boldsymbol{u} = \boldsymbol{Nd_{e}}$

ここで $\boldsymbol{d_{e}}$ は要素の各節点の変位をまとめた節点変位ベクトルであり、各ノードの変位( $(x, y)$ 座標系)を $(u_{i}, v_{i})$ とすると以下です。

$\boldsymbol{d_{e}} = \begin{bmatrix} u_{1} \\ v_{1} \\ u_{2} \\ v_{2} \\ u_{3} \\ v_{3} \\ u_{4} \\ v_{4} \\ \end{bmatrix}$

このように座標と変位に同じ形状関数を用いれる要素をアイソパラメトリック要素といい、計算が非常に簡単になることからFEMではこの要素を主に使用します。

仮想仕事の原理の式に代入

それでは最初に挙げた3つの式に座標変換の式を代入してみましょう。

まずはひずみと変位の関係式です。

$\boldsymbol{\epsilon} = \boldsymbol{LNd_{e}} = \boldsymbol{Bd_{e}}$

ここでBマトリクスと呼ばれる $\boldsymbol{B}$ は以下の形になります。

$\boldsymbol{B} = \begin{bmatrix} \dfrac{\partial N_{1}}{\partial x} & 0 & \dfrac{\partial N_{2}}{\partial x} & 0 & \dfrac{\partial N_{3}}{\partial x} & 0 & \dfrac{\partial N_{4}}{\partial x} & 0 \\ 0 & \dfrac{\partial N_{1}}{\partial y} & 0 & \dfrac{\partial N_{2}}{\partial y} & 0 & \dfrac{\partial N_{3}}{\partial y} & 0 & \dfrac{\partial N_{4}}{\partial y} \\ \dfrac{\partial N_{1}}{\partial x} & \dfrac{\partial N_{1}}{\partial y} & \dfrac{\partial N_{2}}{\partial x} & \dfrac{\partial N_{2}}{\partial y} & \dfrac{\partial N_{3}}{\partial x} & \dfrac{\partial N_{3}}{\partial y} & \dfrac{\partial N_{4}}{\partial x} & \dfrac{\partial N_{4}}{\partial y} \\ \end{bmatrix}$

$N$ は $(\xi, \eta)$ の関数で表されるため、 $(x, y)$ で偏微分できません。ヤコビアン $\boldsymbol{J}$ を使って以下のように変形して求めます。

$\begin{bmatrix} \dfrac{\partial N_{i}}{\partial x} \\ \dfrac{\partial N_{i}}{\partial y} \\ \end{bmatrix} = \begin{bmatrix} \dfrac{\partial x}{\partial \xi} & \dfrac{\partial y}{\partial \xi}\\ \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \eta}\\ \end{bmatrix}^{-1} \begin{bmatrix} \dfrac{\partial N_{i}}{\partial \xi} \\ \dfrac{\partial N_{i}}{\partial \eta} \\ \end{bmatrix} = \boldsymbol{J}^{-1} \begin{bmatrix} \dfrac{\partial N_{i}}{\partial \xi} \\ \dfrac{\partial N_{i}}{\partial \eta} \\ \end{bmatrix}$

次は応力とひずみの関係式です。

$\boldsymbol{\sigma} = \boldsymbol{DBd_{e}}$

最後に仮想仕事の原理の式です。

$\boldsymbol{f_{e}} = h \iint \boldsymbol{B}^{T}\boldsymbol{DB}dxdy \boldsymbol{d_{e}}$

二重積分の部分は座標変換したのちにガウスルジャンドル積分を用いて以下のように書き表せます。
この辺りは実際に要素剛性マトリクスを求めるときに詳しくやります。

${\displaystyle \boldsymbol{f_{e}} = h \sum_{i=1}^{2} \sum_{j=1}^{2} w_{i}w_{j}\boldsymbol{B}^{T}(\xi_{i}, \eta_{i})\boldsymbol{DB}(\xi_{i}, \eta_{i})\mathrm{det}\boldsymbol{J}d\xi d\eta \, \boldsymbol{d_{e}} }$