平面トラスの有限要素法のプログラムをPythonで組んでみました。以前に自分がつくったコードを引っ張り出してきて、Pythonに移植しました。このくらいのボリュームのコードにじっくり取り組むと面白いですね。
トラスのプログラムって実際に使うことはあまりないですが、最初にFEMのコードを組むときはやはり平面トラスからですね。平面トラスが組めれば、有限要素法の計算の流れがわかるので、立体フレームなどへの改良は多少式が複雑になるだけです。
１．サンプルコードについて
以下、サンプルコードの書式について概要を説明したいと思います。
(1) 可読性
コードの可読性を高めるために、def関数へのデータの流れをできるだけ分かりやすくする必要があります。長めのコードは組んだ直後には読めても、１ヶ月ぐらいすると「なんだっけ？」となります。トラスのコードをベースにフレームのコードを組みたいので、ここでの手抜きはあとにひびきます。
さて、有限要素法の古い本はFORTRANで書かれたものが多いので、図書館でFORTRANのコードと格闘することがあります。FORTRANは文法がシンプルなので気長に読めばなんとかなりますが、まるで古文書を解読しているような気分になれます。Pythonとの文法的な違いはFORTRANは引数で渡さなければ基本的に関数の外の変数は参照できません。また、引数は参照渡しなのですが、関数内での引数の内容変更は例外なく関数の外にも影響を与えます。つまり、データの引き渡しは特殊な宣言(COMMON, GLOBAL)をしなければ、入力も出力も基本的には引数で行います。
プログラム言語によってデータの流れが微妙に違うので、移植する場合にはどのようにすれば可読性が高くなるのか試行錯誤することになります。
Pythonでは、関数へのデータの引き渡しにおいて外側の変数を引数なしで内側から参照できるため、関数の引数をほとんどなしにすることも可能です。でも、データの流れはわかりにくくなるように思います。
次に、Pythonの引数では参照値を渡しますが、引数で渡した変数の変更が外側に影響を与える場合と与えない場合があるので注意が必要です。なお、データの受け渡しにglobal変数を使えば楽ですが、不用意に多用すると読みにくいコードになってしまいます。
そこで、私は以下のような方針でコードを記述することにしました。
- 入力データを扱うdef datum()、計算の準備データを扱うbc_index()では、あまりにも関数から読み出すデータが多いため、関数内のデータはすべてglobal変数として関数の外に渡す。
- 関数に渡すデータは、文法上不要であっても引数の部分に記述する。例えばサンプルコードのdef k_elem(iii, indv, xx, eyg, are)では、全ての引数は文法上不要ですが、可読性を高めるためだけに記述する。
- 関数の計算結果はreturanで返す。Pythonでは2つ以上の変数でもreturanで返すことができるので、かなり柔軟に対応できます。
- global変数は不用意には使わない。
上記の方針から外れたものは全体剛性マトリクスのgk、荷重ベクトルf_loadです。特にgkは大きな配列となるため、関数を呼び出す度にローカル変数としてデータが新規に定義されると計算速度に影響するためglobal変数としました。f_loadはglobal変数ではなくてもよかったのかな・・・
(2) その他の書式についての備考
- あとで読み返すときに時間がかかるのが嫌いなのでコメント文は多めです。多分、多すぎると思います。本来は解説資料を別に作って、プログラムはシンプルにするのがよいのだと思います。
- 実用的な使用方法として、ファイルで入出力を行うことになると思いますが、トラスプログラムはFEMの勉強が主な目的だと思いますので、小回りがきくように関数にデータを直接入力するようにしています。
- 出力が美しくありませんが、少数（float）の書式で桁数の指定は行っていません。これは、コードを試行錯誤する上で小数点以下の細かい誤差がわかるようにしたかったからです。
- 最も解析に時間がかかるのは連立一次方程式の部分です。このサンプルコードではnumpy.linalg.solve()を使っています。numpyは計算速度に定評のあるライブラリであり、たった１行で解が求まるのは本当に楽です。なお、規模の大きい構造物など解析時間が気になる場合には、反復法や、疎行列であることを利用した計算方法を考えてみるとよいと思います。
- 入力データ例のdef datum()、def datum_exp01()において、ヤング係数eygは計算内容の検証時に数値が単純になるように1.0を入力しています。
２．さっそく、サンプルコード
さっそく、サンプルコードの紹介です。なお、コードの解説は、別の記事で追加していきたいと思います。
# -*- coding: utf-8 -*-
'''平面トラス解析プログラム Ver.1.01 (2019.11.19) By Snow Tree in June
平面トラスの節点変位、部材応力、反力の算定
'''
import numpy as np
def datum():
'''データ入力
[global変数として入力データを受け渡し]
ndim(=2):ディメンジョン。自由度はx方向、y方向の2種類
uint_max(=65535):'uint16' 0～65535
nod:総節点数
nel:総部材数
nfix:境界条件を与える節点の数
nload:外力が作用する節点の数
eyg :ヤング係数
xx[i][j]:i節点の座標(j=0:x方向, j=1:y方向)
are[i] :i部材の断面積
indv[i][j]:i部材を構成する節点番号(j=0, 1)
npfix[i]:境界条件を与えられている節点番号
ndfix[i][j]:fixの方向:0, freeの方向:1(j=0:x方向, J=1:y方向)
npload[i]:外力が作用する節点番号
aload[i][j]:外力の大きさ(j=0:x方向, J=1:y方向)
'''
global ndim, uint_max
global nod, nel, nfix, nload
global eyg
global xx, are, indv, npfix, ndfix, npload, aload
ndim = 2
uint_max = 65535
nod = 3
nel = 3
nfix = 2
nload = 1
eyg = 1.
#'uint16' 0~65535
#numpy配列の宣言
xx = np.zeros((nod, ndim), dtype = 'float64')
are = np.zeros(nel, dtype = 'float64')
indv = np.zeros((nel, ndim), dtype = 'uint16')
npfix = np.zeros(nfix, dtype = 'uint16')
ndfix = np.zeros((nfix, ndim), dtype = 'bool')
npload = np.zeros(nload, dtype = 'uint16')
aload = np.zeros((nload, ndim), dtype = 'float64')
xx[:] = [[-100, 0],
[0, 100],
[100, 0]]
are[:] = 100.
indv[:] = [[0,1],
[1,2],
[0,2]]
npfix[:] = [2,0]
ndfix[:] = [[1, 0],
[0, 0]]
npload[0] = 1
aload[0][:] = [0., -3.]
def datum_exp01():
'''データ入力例
5節点、7部材トラスの例
'''
global ndim, uint_max
global nod, nel, nfix, nload
global eyg
global xx, are, indv, npfix, ndfix, npload, aload
ndim = 2
uint_max = 65535
nod = 5
nel = 7
nfix = 2
nload = 2
eyg = 1.
#'uint16' 0~65535
#numpy配列の宣言
xx = np.zeros((nod, ndim), dtype = 'float64')
are = np.zeros(nel, dtype = 'float64')
indv = np.zeros((nel, ndim), dtype = 'uint16')
npfix = np.zeros(nfix, dtype = 'uint16')
ndfix = np.zeros((nfix, ndim), dtype = 'bool')
npload = np.zeros(nload, dtype = 'uint16')
aload = np.zeros((nload, ndim), dtype = 'float64')
xx[:] = [[0., 0.],
[100., 100.],
[200., 0.],
[300., 100.],
[400., 0.]]
are[:] = 100.
indv[:] = [[0, 1],
[1, 2],
[2, 3],
[3, 4],
[0, 2],
[2, 4],
[1, 3]]
npfix = [0, 4]
ndfix[:] = [[0, 0],
[1, 0]]
npload[:] = [1, 3]
aload[:] = [[0., -3.],
[0., -3.]]
def output_datum():
'''入力データの出力
'''
print('========================================')
print('１．入力データ')
print('(1) 制御データ')
print(' ndim, uint_max：[{}, {}]'.format(ndim,uint_max))
print('(2) 共通データ')
print(' E：[{}]'.format(eyg))
print('(3) 節点座標')
print(' 節点数：{}'.format(nod))
for i in range(nod):
print('{:4d}: [x:{:12.5f}, y:{:12.5f}]'\
.format(i, xx[i][0], xx[i][1]))
print('(4) 部材')
print(' 部材数：{}'.format(nel))
for i in range(nel):
print('{:4d}: (node:{:4d} -- {:4d}), [A:{:12.5f}]'\
.format(i, indv[i][0], indv[i][1], are[i]))
print('(5) 境界条件')
print(' 境界条件を与える節点の数：{}'.format(nfix))
print(' [境界条件：fix:0, free:1]')
for i in range(nfix):
print('{:4d}: (node:{:4d}), [{:2d}, {:2d}]'\
.format(i, npfix[i], ndfix[i][0], ndfix[i][1]))
print('(6) 外力')
print(' 外力が作用する節点の数：{}'.format(nload))
for i in range(nload):
print('{:4d}: (node:{:4d}), [x:{:12.5f}, y:{:12.5f}]'\
.format(i, npload[i], aload[i][0], aload[i][1]))
def bc_index(ndim, uint_max, nod, nfix, npfix, ndfix):
'''計算準備データの生成
[global変数として結果を受け渡し]
nmat:節点数(nod) * ディメンジョン(ndim)
n_free:境界条件がfreeの部分の配列の大きさ
n_fix:境界条件がfixの部分の配列の大きさ
ind_free[nmat → n_free]:free部分の通し番号と節点番号の対応表
ind_fix[n_fix]:fix部分の通し番号と節点番号の対応表
'''
global nmat, n_free, n_fix, ind_free, ind_fix
nmat = nod * ndim
if nmat >= uint_max:
print('節点数が制限値を超えている')
n_fix = 0
for i in range(nfix):
for j in range(ndim):
if ndfix[i][j] == 0:
n_fix += 1
n_free = nmat - n_fix
ind_free = np.arange(nmat, dtype = 'uint16')
ind_fix = np.zeros(n_fix, dtype = 'uint16')
#fixの部分をuint_maxに入れ替え
for i in range(nfix):
ii = npfix[i]
for j in range(ndim):
if ndfix[i][j] == 0:
ind_free[ndim * ii + j] = uint_max
#fixの部分を除き、前につめる。[0～mm]の範囲が必要なデータ
s0 = 0
s1 = 0
for i in range(nmat):
if ind_free[i] == uint_max:
ind_fix[s0] = i
s0 += 1
else:
ind_free[s1] = ind_free[i]
s1 += 1
#本来は以下の1行は不要。ind_freeの内容を読みやすくするために削除
ind_free = np.delete(ind_free, slice(n_free, nmat), 0)
def k_elem(iii, indv, xx, eyg, are):
'''iii番目の要素の全体座標系要素剛性マトリクスの生成
[return]
ek_global[4][4]:全体座標系要素剛性マトリクス
'''
ek_local = np.zeros((2, 2), dtype = 'float64')
ek_global = np.zeros((4, 4), dtype = 'float64')
xl, rotation = transform(iii, indv, xx)
#部材座標系要素剛性マトリクス
ek_local[0][0] = eyg * are[iii] / xl
ek_local[0][1] = - ek_local[0][0]
ek_local[1][0] = ek_local[0][1]
ek_local[1][1] = ek_local[0][0]
#全体座標系要素剛性マトリクス
ek_global = np.dot(rotation.T, np.dot(ek_local, rotation))
return ek_global
def transform(iii, indv, xx):
'''indv[iii]の部材長、[[c,s,0,0],[0,0,c,s]]の算定
[return]
xl: 部材長
rotation:[[cos(θ), sin(θ), 0, 0],
[0, 0, cos(θ), sin(θ)]]
（全体座標系における部材の傾きをθとする。反時計回りを正。）
'''
rotation = np.zeros((2,4), dtype ='float64')
k0 = indv[iii][0]
k1 = indv[iii][1]
x21 = xx[k1][0] - xx[k0][0]
y21 = xx[k1][1] - xx[k0][1]
xl = np.sqrt(x21 * x21 + y21 * y21)
c = x21 / xl
s = y21 / xl
rotation[:] = [[c, s, 0., 0.],
[0., 0., c, s]]
return xl, rotation
def set_stiff(iii, ek_global, indv):
'''要素剛性マトリクスから全体剛性マトリクスの生成（fix,freeを含む）
[global変数]
gk[nmat][nmat]:全体剛性マトリクス（fix,freeを含む）
'''
global gk
m = 2
for i in range(m):
i_t = ndim * indv[iii][i]
i_e = ndim * i
for j in range(m):
j_t = ndim * indv[iii][j]
j_e = ndim * j
for ij in range(ndim):
ij_t = i_t + ij
ij_e = i_e + ij
for ji in range(ndim):
ji_t = j_t + ji
ji_e = j_e + ji
gk[ij_t][ji_t] += ek_global[ij_e][ji_e]
return gk
def set_load(ndim, nmat, nload, npload, aload):
'''荷重ベクトルの生成（fix,freeを含む）
[return]
f_load[nmat]:荷重ベクトル
'''
f_load = np.zeros(nmat, dtype = 'float64')
if nload != 0:
for i in range(nload):
k = ndim * npload[i]
for j in range(ndim):
f_load[k + j] += aload[i][j]
return f_load
def extract_fix(nmat, n_free, n_fix, ind_free, ind_fix, gk):
'''全体剛性マトリクスgkのfix部分をk_yfixに代入（反力算定用）
[return]
k_yfix[n_fix][n_free]:全体剛性マトリクスgk[nmat][nmat]のfix部分
'''
k_yfix = np.zeros((n_fix, n_free), dtype = 'float64')
for i in range(n_fix):
ii = ind_fix[i]
for j in range(n_free):
jj = ind_free[j]
k_yfix[i][j] = gk[ii][jj]
return k_yfix
def del_stiff_fix(nmat, n_free, ind_free):
'''f_load[], gk[][]のfix部分を削除
[global変数]
f_load[nmat → n_free]:荷重ベクトル(fix部削除)
gk[nmat → n_free][nmat → n_free]:
全体座標系剛性マトリクス(fix部削除)
'''
global f_load, gk
for i in range(n_free):
ii = ind_free[i]
f_load[i] = f_load[ii]
for j in range(n_free):
jj = ind_free[j]
gk[i][j] = gk[ii][jj]
#行れ鵜計算をするために削除。※numpy配列の形状変更は計算時間を要する
f_load = np.delete(f_load, slice(n_free, nmat), 0)
gk = np.delete(gk, slice(n_free, nmat), 0)
gk = np.delete(gk, slice(n_free, nmat), 1)
def calc_disp(ndim, nod, n_free, ind_free, y_free):
'''節点変位の計算
[return]
disp_c[i][j]:i節点のj方向の変位
(j=0:x方向変位, J=1:y方向変位, j=2:変位の大きさ)
'''
disp_c = np.zeros((nod, ndim + 1), dtype = 'float64')
for i in range(n_free):
ii = ind_free[i] // ndim
jj = ind_free[i] % ndim
disp_c[ii][jj] = y_free[i]
for i in range(nod):
dis_x = disp_c[i][0]
dis_y = disp_c[i][1]
disp_c[i][2] = np.sqrt(dis_x ** 2 + dis_y ** 2)
return disp_c
def calc_stress(nel, ndim, indv, xx, disp_c):
'''部材軸力の計算
[return]
stress_c[i]:i部材の軸力（引張側正）
'''
em_tmp = np.zeros(4, dtype = 'float64')
stress_c = np.zeros(nel, dtype = 'float64')
for iii in range(nel):
xl, rotation = transform(iii, indv, xx)
for i in range(2):
k1 = indv[iii][i]
for j in range(ndim):
k2 = ndim * i + j
em_tmp[k2] = disp_c[k1][j]
xs = np.dot(rotation, em_tmp)
stress_c[iii] = eyg * are[iii] * (xs[1] - xs[0]) / xl
return stress_c
def calc_react_f(n_fix, ind_fix, k_yfix, y_free):
'''反力の計算
[return]
rf_c[i][j]:(0 ≦ i ＜ n_fix)
(j=0:反力の節点番号, j=1:反力の向き(x:0, y:1), j=2:反力)
'''
rf_c = np.dot(k_yfix, y_free)
rf_dere = [[],[]]
for i in range(n_fix):
rf_dere[0].append(ind_fix[i] // ndim)
if (ind_fix[i] % ndim) == 0:
rf_dere[1].append('x')
else:
rf_dere[1].append('y')
return rf_c, rf_dere
def output_result():
'''計算結果の出力
誤差確認のためfloatの出力は書式指定せず
'''
print('========================================')
print('２．計算結果')
print('(1) 節点変位')
print(' 節点: 変位[x:], [y:], [xy:]')
for i in range(nod):
print('{:4d}: [x:{}], [y:{}], [xy:{}]'\
.format(i, disp_c[i][0], disp_c[i][1], disp_c[i][2]))
print('(2) 部材応力')
print(' 部材: (節点1--節点2), 部材応力')
for i in range(nel):
print('{:4d}: ({:4d} -- {:4d}), {}'\
.format(i, indv[i][0], indv[i][1], stress_c[i]))
print('(3) 反力')
print(' No.:(節点-方向), 反力')
for i in range(n_fix):
print('{:4d}:({:4d} - {}), {}'\
.format(i, rf_dere[0][i], rf_dere[1][i], rf_c[i]))
if __name__ == '__main__':
'''メインプログラム
[入力、計算準備データ]
datum(), bc_index() にて
[主な計算用データ]
gk[nmat→nfree]:全体剛性マトリクス
f_load[nmat→nfree]:荷重ベクトル
ek_local[2][2]:要素剛性マトリクス（部材座標系）
ek_global[4][4]：要素剛性マトリクス（全体座標系）
k_yfix[n_fix][n_free]:変位がfixされた部分の剛性マトリクスの抜き出し
y_free[n_free]：節点変位のベクトル（freeの部分）
[計算結果]
disp_c[nod]:節点変位のベクトル（free,fix）
stress_c[nel]：部材応力
rf_c[n_fix]：反力
rf_dere[n_fix][i]：反力の節点番号と方向(i=0:節点番号、i=1:方向)
'''
#データ入力
datum()
#datum_exp01()
#入力データの出力
output_datum()
#計算準備データ
#global nmat, n_free, n_fix, ind_free, ind_fix
bc_index(ndim, uint_max, nod, nfix, npfix, ndfix)
gk = np.zeros((nmat, nmat), dtype = 'float64')
for iii in range(nel):
#全体座標系要素剛性マトリクスの生成
ek_global = k_elem(iii, indv, xx, eyg, are)
#要素剛性マトリクスより全体剛性マトリクスを生成
#global gk
set_stiff(iii, ek_global, indv)
#荷重ベクトルの生成
f_load = set_load(ndim, nmat, nload, npload, aload)
#fixされているところの剛性マトリクスのみを抽出（反力算定用）
k_yfix = extract_fix(nmat, n_free, n_fix, ind_free, ind_fix, gk)
#f_load[], gk[][]のfixの部分を削除
#global f_load, gk
del_stiff_fix(nmat, n_free, ind_free)
#連立一次方程式を解き、拘束条件のない部分の変位を求める
y_free = np.linalg.solve(gk, f_load)
#節点変位の計算
disp_c = calc_disp(ndim, nod, n_free, ind_free, y_free)
#部材応力の計算
stress_c = calc_stress(nel, ndim, indv, xx, disp_c)
#反力の計算
rf_c, rf_dere =calc_react_f(n_fix, ind_fix, k_yfix, y_free)
#計算結果の出力
output_result()私が実際に購入した教材のご紹介
以下、私が実際に購入したPythonの教材をまとめてみました。 Pythonを学習する上で、少しでもお役に立つことができればうれしいです。
