Pythonによる2次元FEM応力解析プログラム


プログラムとしての仕上がりは下の投稿のほうがいいと思いますので参考にしていただければと思います．

+ [Pythonによる２次元応力解析プログラム（四角形要素）の理論を少し詳しく解説する（改訂版）](https://qiita.com/damyarou/items/320bad2052bb5fccd75f)
+ [有限要素法とPythonプログラム（三角形一次要素を用いた２次元応力解析プログラム）](https://qiita.com/damyarou/items/8ca3432f9bba20a1f5ea)

## 概要
Pythonによる線形弾性２次元FEM応力解析プログラムを作ってみた．

要素は，１要素４節点のアイソパラメトリック要素であり，１要素に４ガウス積分点を持つ．

扱える荷重は以下の通り．

+ 節点外力
+ 節点温度変化
+ 自重および地震時慣性力

### 理論式
##### 節点変位計算式
$$
[\boldsymbol{k}]\\{\boldsymbol{u}\\}=\\{\boldsymbol{f}\\}+\\{\boldsymbol{f_t}\\}+\\{\boldsymbol{f_b}\\}
$$

$$
\begin{align}
&[\boldsymbol{k}]=t\cdot \int_A [\boldsymbol{B}]^T[\boldsymbol{D}][\boldsymbol{B}] dA \\\
&\\{\boldsymbol{f_t}\\}=t\cdot \int_A [\boldsymbol{B}]^T[\boldsymbol{D}]\\{\boldsymbol{\epsilon_0}\\} dA \\\
&\\{\boldsymbol{f_b}\\}=t\cdot\gamma\cdot \int_A [\boldsymbol{N}]^T [\boldsymbol{N}] dA \cdot \\{\boldsymbol{w}\\}
\end{align}
$$

+ $[\boldsymbol{k}]$ : 要素剛性マトリックス
+ $\\{\boldsymbol{u}\\}$ : 節点変位ベクトル
+ $\\{\boldsymbol{f}\\}$ : 節点外力ベクトル
+ $\\{\boldsymbol{f_t}\\}$ : 温度変化による節点力ベクトル
+ $\\{\boldsymbol{f_b}\\}$ : 自重・地震加速度による節点ベクトル
+ $[\boldsymbol{D}]$ : 応力ー歪関係マトリックス
+ $t$, $\gamma$ : 要素厚さ，要素単位体積重量
+ $\\{\boldsymbol{w}\\}$ : 節点加速度ベクトル（重力加速度に対する比率で入力）

##### 要素応力算定式
$$
\\{\boldsymbol{\sigma}\\}=[\boldsymbol{D}]\\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\\}
$$
$\\{\boldsymbol{\epsilon_0}\\}$は，温度変化による歪

##### 節点変位ー任意位置ひずみ関係式
$$
\\{\boldsymbol{\epsilon}\\}=[\boldsymbol{B}]\\{\boldsymbol{u}\\}
$$
$\\{\boldsymbol{\epsilon}\\}$は，要素内任意点の歪

##### 節点変位ー任意位置歪算定式
$$
\\{\boldsymbol{v}\\}=[\boldsymbol{N}]\\{\boldsymbol{u}\\}
$$
$\\{\boldsymbol{v}\\}$は，要素内任意点の変位

### 使えるか？
Pythonはインタプリタなので，計算速度的に使い物になるか不安であったが，節点数：3417，要素数：3257，自由度：6834のモデルでの計算時間は，７．９秒であり，使えるレベルであることが確認できた．

## プログラム
プログラムは，長いので，Gistに貼り付けたものにリンクを貼っています．
なお，コンター作図といっても，要素を応力レベルに応じて色分けして塗りつぶしているだけです．私の目的では別に内挿してコンターラインを引いて貰う必要はないため．

+ [２次元FEM応力解析プログラム（py_fem_plnt.py）](https://gist.github.com/damyarou/9ed70fd1365177fe1b6fcf49f7b1aecc/raw/1_py_fem_plnt.py)
+ [主応力コンター作図プログラム（py_fig_cont.py）](https://gist.github.com/damyarou/9ed70fd1365177fe1b6fcf49f7b1aecc/raw/2_py_fig_cont.py)
+ [主応力ベクトル＆変位モード作図プログラム（py_fig_vect.py）](https://gist.github.com/damyarou/9ed70fd1365177fe1b6fcf49f7b1aecc/raw/3_py_fig_vect.py)
+ [3417節点モデルの入力データ（inp_arch.txt）](https://gist.github.com/damyarou/9ed70fd1365177fe1b6fcf49f7b1aecc/raw/4_inp_arch.txt)


## 入力データフォーマット

```
npoin  nele  nsec  npfix  nlod  NSTR # 基本量
t  E  po  alpha  gamma  gkh  gkv     # 材料特性
..... (1~nsec) .....
node1  node2  node3  node4  isec     # 要素-節点関係，材料特性番号
..... (1~nele) .....
x  y  deltaT                         # 節点座標，節点温度変化
..... (1~npoin) .....
node  kox  koy  rdisx  rdisy         # 変位拘束条件
..... (1~npfix) .....
node  fx  fy                         # 外力
..... (1~nlod) .....
```

<table>
<tr><td>npoin, nele, nsec</td><td>節点数，要素数，材料特性数</td></tr>
<tr><td>npfix, nlod</td><td>拘束節点数，載荷節点数</td></tr>
<tr><td>NSTR</td><td>応力状態（平面歪：０，平面応力：１）</td></tr>
<tr><td>t, E, po, alpha</td><td>板厚，弾性係数，ポアソン比，線膨張係数</td></tr>
<tr><td>gamma, gkh, gkv</td><td>単位体積重量，水平及び鉛直方向加速度（gの比）</td></tr>
<tr><td>x, y, deltaT</td><td>節点x座標，節点y座標，節点温度変化</td></tr>
<tr><td>node, kox, koy</td><td>拘束節点番号，x及びy方向拘束の有無（拘束：１，自由：０）</td></tr>
<tr><td>rdisx, rdisy</td><td>x及びy方向変位（無拘束でも０を入力）</td></tr>
<tr><td>node, fx, fy</td><td>載荷重節点番号，x方向荷重，y方向荷重</td></tr>
</table>

## 出力データフォーマット

```
npoin  nele  nsec npfix  nlod   NSTR
    4     1     1     2     2      1
  sec               t               E              po           alpha           gamma        gkh        gkv
    1   1.0000000e+00   1.0000000e+03   0.0000000e+00   1.0000000e-05   2.3000000e+00      0.000      0.000
 node               x               y              fx              fy          deltaT   kox   koy
    1   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     1     1
    2   1.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00     0     1
    3   1.0000000e+00   1.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     0     0
    4   0.0000000e+00   1.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     0     0
 node   kox   koy          rdis_x          rdis_y
    1     1     1   0.0000000e+00   0.0000000e+00
    2     0     1   0.0000000e+00   0.0000000e+00
 elem     i     j     k     l   sec
    1     1     2     3     4     1
 node           dis-x           dis-y
    1   0.0000000e+00   0.0000000e+00
    2   8.8817842e-19   0.0000000e+00
    3   1.7763568e-18   2.0000000e-02
    4   4.1448326e-18   2.0000000e-02
 elem           sig_x           sig_y          tau_xy              p1              p2             ang
    1  -7.4014868e-16   2.0000000e+01   1.1333527e-16   2.0000000e+01   0.0000000e+00   9.0000000e+01
n=8  time=0.044 sec
```

<table>
<tr><td>node, dis-x, dis-y</td><td>節点番号，x方向変位，y方向変位</td></tr>
<tr><td>elme, sig_x, sig_y, tau_xy</td><td>要素番号，x方向直応力，y方向直応力，せん断応力</td></tr>
<tr><td>p1, p2, ang</td><td>第一主応力，第二主応力，第一主応力の方向</td></tr>
<tr><td>n, time</td><td>全自由度，計算時間</td></tr>
</table>

## 出力事例

![tex_fig.png](https://qiita-image-store.s3.amazonaws.com/0/129300/bb93e101-0728-7b24-dfc0-c9c36a2fb08a.png)



以　上