マトリクス構造解析オンライン授業用テキスト
第12回オンライン授業
YouTube授業動画リスト
小さい字は補足説明なので、読み飛ばしてもいいです。
今回は、プログラムを使ってトラスの問題を解いてみよう。 第5回で2要素について手計算で解いたトラスの 問題を、もっと部材数の多い普通のトラスについて、 プログラムで解いてみようという趣旨だ。 本来であれば、プログラムを組むところからやってもらうのが望ましいが、 この授業の枠内で、これぐらいのプログラムを組めるようになってもらうのは、 なかなか難しいので、ここでは、こちらが用意したプログラムを使って 解いてもらうことにする。
では、このtorasu2.pyの使い方を説明するために、
第5回の2要素トラスを
torasu2.pyで解く手順を説明する。
節点変位の正解をメモしておく。
$v_{2}=(1+2\sqrt{2})\frac{P\ell}{EA}=(1+2\sqrt{2})\frac{P}{k}$
$w_{2}=\frac{P\ell}{EA}=\frac{P}{k}$
torasu2.pyをテキストエディター(geditなど)で開き、 上の方の、以下のように書いてあるところを見つける。
nyou = 2 #要素数
nset = 3 #節点数
nyouは要素数、nsetは節点数を表している。 上記の2要素トラスの例題は、要素数が2要素で、 節点数が3個だから、nyou=2 で、nset=3 となる。 他の要素数、節点数の問題を解く場合は、ここを書き換える。
次に、もう少し下の以下のように書いてあるところを見つける。
E = 1.0 #ヤング率 N/mm^2(MPa)
ell = 1.0 #1要素の長さ mm
A = 1.0 #面積 mm^2
sk1 = E*A/ell #要素1のばね定数
# 要素剛性マトリクス
s[1][2][2] = sk1; s[1][2][4] = -sk1
s[1][4][2] = -sk1; s[1][4][4] = sk1
sk2 = sk1 / math.sqrt(2.0) #要素2のばね定数
s[2][2][2] = sk2; s[2][2][4] = -sk2
s[2][4][2] = -sk2; s[2][4][4] = sk2
コマンドや式を1行に並べて書きたいときは ; で区切る。
Eはヤング率、ellは要素の長さ、aは断面積、sk1は、
要素①の ばね定数$k_{①}$を表している。
プログラムで各要素の要素剛性方程式は、
その下に書かれているs[1][2][2]みたいな配列で与えている。
s[1][2][2]の最初の[1]が、要素①という意味で、
次の[2][2]は、$k_{①}$の2行2列の成分$k^{①}_{2,2}$を表している。
ここで、各要素の局所系の要素剛性方程式を与えていきたいのだが、
各要素の局所系の要素剛性方程式は、
第4回でやったように
以下のような形をしている。
$
\left(
\begin{array}{c}
S_{1}^{\ell}\\
N_{1}^{\ell}\\
S_{2}^{\ell}\\
N_{2}^{\ell}
\end{array}
\right)
=
\left(
\begin{array}{rrrr}
0 & 0 & 0 & 0 \\
0 & k & 0 & -k \\
0 & 0 & 0 & 0 \\
0 & -k & 0 & k \\
\end{array}
\right)
\left(
\begin{array}{c}
v_{1}^{\ell}\\
w_{1}^{\ell}\\
v_{2}^{\ell}\\
w_{2}^{\ell}
\end{array}
\right)
$
つまり、各要素について、2行2列と4行4列に、その要素のばね定数$k$を、
2行4列と、4行2列に、$-k$を与える。
その他の成分は、あらかじめ初期化しているので、
改めて0を与えなくてもよい。
例えば、要素①の ばね定数を$k_{①}=\frac{EA}{\ell}$で与えるとすると、
$E$や$A$や$\ell$を
E=10.e3[N/mm^2], ell=1000.0[mm], a=10.0**2[mm^2]みたいに与えておいて、
それからsk1=E*a/ellとすれば、sk1に$\frac{EA}{\ell}$が代入されるから、
これを要素①の剛性方程式の$k_{①}$として、
s[1][2][2] = sk1; s[1][2][4] = -sk1
みたいに与えていけばよい。
さて、プログラムで計算するには、$E$や$A$や$\ell$には具体的な値を 入れなければならない。 今回は、プログラムがちゃんと動いているかどうかの確認のため、 $E=A=1$, 要素①について$\ell_{①}=1$, 要素②について$\ell_{②}=\sqrt{2}$として計算してみる。 プログラムの動作確認の際に、 各種のパラメータをとりあえず全部「1」にして計算してみるという ことはよくやるが、これはこれで注意が必要である。 ヤング率も断面積も長さも全部「1」という構造は、 実際には、ものすごく柔らかいフニャフニャの材料だったり、 逆にものすごく固い剛体のような材料だったり、 ありえない構造を解いていることになる。 梁の弾性計算程度であれば、そういう変な値で解いても、 まあ、正解が求まるかもしれないが、 あまりに極端な材料であるために、計算できなくなったり、 おかしな値が出ることもある。 そういうことを避けるためには、 自分が想像しやすい大きさの現実的な材料で解いてみるのがよい。 例えば、 E=10.e3[N/mm^2] は木材程度のヤング率で、 10mm$\times$10mmの角材ぐらいの太さ(a=10.0**2[mm^2])で、 1m(ell=1000.0[mm])ぐらいの部材を想定してみれば、 10mm角の長さ1mの角材でトラスを組んで、 人間1人(50kgfとか)ぐらいがぶらさがったら、何mmぐらいたわみそうかとか、 ある程度の予測ができる。 そういう予測のできる問題で解いてみるということも重要だ。 なのだが、今回は、出てきた数値が単位荷重法の理論値と比較しやすいことを 優先して $E=A=1$, 要素①について$\ell_{①}=1$, 要素②について$\ell_{②}=\sqrt{2}$として計算してみる。 つまり、$k_{①}=\frac{EA}{\ell}=1$, $\;\;k_{②}=\frac{1}{\sqrt{2}}$として 計算する ということである。 まず、要素①について、以下のように E=1.0, ell=1.0, a=1.0と与えてやれば、 sk1=E*a/ellは1になる。
E = 1.0 #ヤング率 N/mm^2(MPa)
ell = 1.0 #1要素の長さ mm
A = 1.0 #面積 mm^2
sk1 = E*A/ell #要素1のばね定数
# 要素剛性マトリクス
s[1][2][2] = sk1; s[1][2][4] = -sk1
s[1][4][2] = -sk1; s[1][4][4] = sk1
次に、その下の要素2の剛性マトリクスのところで、 $k_{②}$を与えているsk2には、以下のように$\frac{k_{①}}{\sqrt{2}}$を与える。 math.sqrt(2.0)は、$\sqrt{2.0}$のことである。
sk2 = sk1 / math.sqrt(2.0) #要素2のばね定数
s[2][2][2] = sk2; s[2][2][4] = -sk2
s[2][4][2] = -sk2; s[2][4][4] = sk2
今回は2要素しかないので、要素3以降の剛性マトリクスは 与える必要がない。
次に各要素の節点番号を与える。 剛性マトリクスの下の以下のように書いているところを見つける。
#各要素の左節点番号、右節点番号
hidari[1] = 1; migi[1] = 2
hidari[2] = 2; migi[2] = 3
hidari[1]の[1]は、要素①ということで、
hidariというのは、左節点番号ということである。
migi[1]というのは、要素①の右節点番号ということである。
例えば、今回の2要素トラスについて、
まずは、図のようにすべての要素を$z$軸方向に並べて初期状態を決める。
トラスが一筆書きで書けるなら、なるべく1本につなげた方が、
設定しやすい。
要素①の左節点は1で右節点は2だから、
hidari[1]=1; migi[1]=2 となる。
要素②の左節点は2で右節点は3だから、
hidari[2]=2; migi[2]=3（これも、このまま修正不要）。
今回は、2要素なので、要素③以降はない。
次に各要素の回転角を与える。 節点番号の下の以下のように書いてある場所を見つける。
#各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th[1] = 0.0
th[2] = pi + pi/4.0
th[1]というのは、要素①の初期状態からの回転量を、
左節点の$x$軸右ねじ回り($yz$平面では反時計回り)の回転角$\theta_{①}$で与える。
要素①は、水平で初期状態から回転しないので、$\theta_{①}=0$である。
つまり、th[1]=0.d0でよい。
一方、要素②は、斜めになっているので、回転角$\theta_{②}$は0ではない。
要素②の左節点(つまり節点2)の
初期状態からの$x$軸右ねじ回りの回転角は、
図のように$\pi+\frac{\pi}{4}$である。
よって、th[2]=pi+pi/4.0と書き換える必要がある。
次は境界条件である。 回転角の下の以下のように書いてあるところを見つける。
# 境界条件
xv[1] = 0.0; xw[1] = 0.0
xv[3] = 0.0; xw[3] = 0.0
境界条件は、拘束されている節点の拘束されている変位 ($v$や$w$)を0にすることで与える。 $v_{1}$はxv[1], $w_{2}$はxw[2]のように表す。 例えば、節点3で$v_{3}=0$であれば、xv[3]=0と与える。 今回は、節点1と節点3が固定端で、 $v_{1}=w_{1}=v_{3}=w_{3}=0$だから、 上のようになる。
最後に載荷条件である。境界条件のちょっと下の 以下のような場所を見つける。
# 荷重条件(y方向:fy[節点番号], z方向:fz[節点番号])
fy[2] = 1.0 #[N]
荷重は、$y$方向荷重がfy[], $z$方向荷重がfz[]で、 荷重を与える節点番号を[]内に書く。 剛性マトリクスのばね定数を与えるところで述べたように、 荷重も、現実的な荷重を与えるべきだが、 今回は、理論値との比較がしやすいように$P=1$ということにする。 節点2の$y$方向に1だから、fy[2]=1.0とすればよい。
プログラムの修正が終わったら、 "torasu2.py"に上書き保存してもいいが、 エラーが出たときに、エラーが出ないプログラムを少しずつ修正しながら、 どの段階でエラーが出るのか比較できるように、 元のプログラムは"torasu2.py"のまま保管しておいて、 修正したプログラムは名前を変えて、例えば"torasu2_gotou.py"みたいに して保存した方がいいかもしれない。 プログラムを実行すると、 以下のように表示される。
節点番号：1 v=0.0000000000, w=0.0000000000
節点番号：2 v=3.8284271247, w=1.0000000000
節点番号：3 v=0.0000000000, w=0.0000000000
部材番号： 1 S(1)= 0.0000000000, N(1)= -1.0000000000
S(2)= 0.0000000000, N(2)= 1.0000000000
部材番号： 2 S(2)= 1.0000000000, N(2)= -1.0000000000
S(3)= -1.0000000000, N(3)= 1.0000000000
これは、節点ごとの節点変位が並んでいる。
上から順に1行目が$v_{1}, w_{1}$, 2行目が$v_{2}, w_{2}$, 3行目が$v_{3}, w_{3}$である。
境界条件から、$v_{1}=0, w_{1}=0, v_{3}=0, w_{3}=0$はいいだろう。
第5回で単位荷重法で解いた答えによると、
$v_{2}=(1+2\sqrt{2})\frac{P\ell}{EA}$,
$\;w_{2}=\frac{P\ell}{EA}$だったから、$P=E=A=\ell=1$を代入すると、
$v_{2}=(1+2\sqrt{2})=3.828427127124746\;\;\;$,
$w_{2}=1$となる。
$v_{2}$も$w_{2}$も、表示されていいる桁はすべて合っている。。
次に節点力だが、
以下の部材番号: 1は、要素①を、部材番号: 2は要素②を表し、
S(1)やN(1)は$S_{1}$や$N_{1}$を表す。つまり、
部材番号: 1のところのS(2)は$S_{2}^{①}$を表し、
部材番号: 2のところのS(2)は$S_{2}^{②}$を表す。
部材番号： 1 S(1)= 0.0000000000, N(1)= -1.0000000000
S(2)= 0.0000000000, N(2)= 1.0000000000
部材番号： 2 S(2)= 1.0000000000, N(2)= -1.0000000000
S(3)= -1.0000000000, N(3)= 1.0000000000
第5回でマトリクスを手計算で解いた答えによると、
要素①の要素剛性方程式に節点変位を代入
$
\left(
\begin{array}{c}
S_{1}\\
N_{1}\\
S_{2}^{①}\\
N_{2}^{①}
\end{array}
\right)
=
\left(
\begin{array}{rrrr}
0 & 0 & 0 & 0 \\
0 & k & 0 & -k \\
0 & 0 & 0 & 0 \\
0 & -k & 0 & k
\end{array}
\right)
\left(
\begin{array}{c}
0\\
0\\
(1+2\sqrt{2})\frac{P}{k}\\
\frac{P}{k}
\end{array}
\right)
=
\left(
\begin{array}{c}
0\\
-P\\
0 \\
P
\end{array}
\right)
$
要素②の要素剛性方程式に節点変位を代入
$
\left(
\begin{array}{c}
S_{2}^{②}\\
N_{2}^{②}\\
S_{3}\\
N_{3}
\end{array}
\right)
=
\frac{k}{2\sqrt{2}}
\left(
\begin{array}{rrrr}
1 & -1 & -1 & 1 \\
-1 & 1 & 1 & -1 \\
-1 & 1 & 1 & -1 \\
1 & -1 & -1 & 1
\end{array}
\right)
\left(
\begin{array}{c}
(1+2\sqrt{2})\frac{P}{k}\\
\frac{P}{k}\\
0\\
0
\end{array}
\right)
=
\left(
\begin{array}{c}
P\\
-P\\
-P \\
P
\end{array}
\right)
$
と求まっていたが、$P=1$を代入すれば、
プログラムで解いた答えと一致していることが確認できる。
以下の2要素トラス（自分の学籍番号のもの）の 載荷点の節点変位($v, w$)を単位荷重法で求めよ。 部材の伸び剛性はすべて$EA$とする。 斜めの部材は水平方向から$45^{\circ}$傾いていて、 その他の部材は、水平か鉛直である。 次に、$P=E=A=\ell=1.$として、 それをtorasu2.pyで解いて、 手計算の答えとどれくらい合うか考察せよ。 WebClassの「プログラム課題4」から、 単位荷重法の手計算を撮影した画像と プログラム本体（例えば torasu2_gotou.py とか）と、 出力ファイル（例えば kekka.out とか）を それぞれ提出せよ。 どうしてもエラーが取れないなどの場合は、 kekka.outの代わりに、エラー画面等 をテキストエディターに貼り付けて保存したファイルを提出する。 また、最後のコメント欄に、手計算とプログラムの計算とを比較した考察や、 うまく実行できなかった場合は、 どういう症状で実行できないのかの説明を書き込むこと。
学籍番号の末尾が1か6の人。
学籍番号の末尾が2か7の人。
学籍番号の末尾が3か8の人。
学籍番号の末尾が4か9の人。
学籍番号の末尾が5か0の人。
import math
def mxmx(a, b, c):
# c = a @ b
for i in range(1, 5):
for j in range(1, 5):
c[i][j] = 0.0
for k in range(1, 5):
c[i][j] += a[i][k] * b[k][j]
def zahyou(t, theta):
# 4x4座標変換行列
for i in range(1, 3):
for j in range(1, 3):
t[i+2][j] = 0.0
t[i][j+2] = 0.0
t[1][1] = math.cos(theta)
t[1][2] = math.sin(theta)
t[2][1] = -math.sin(theta)
t[2][2] = math.cos(theta)
t[3][3] = t[1][1]
t[3][4] = t[1][2]
t[4][3] = t[2][1]
t[4][4] = t[2][2]
def mxtmx(a, b, c):
# c = a^T @ b
for i in range(1, 5):
for j in range(1, 5):
c[i][j] = 0.0
for k in range(1, 5):
c[i][j] += a[k][i] * b[k][j]
def gausu(n, a, x, b):
# ガウス消去法による連立一次方程式 a x = b の解
# a: (n+1)x(n+1), x: (n+1), b: (n+1) (0番目ダミー)
for k in range(1, n):
for i in range(k+1, n+1):
for j in range(k+1, n+1):
a[i][j] -= a[k][j] * a[i][k] / a[k][k]
b[i] -= b[k] * a[i][k] / a[k][k]
# 後退代入
x[n] = b[n] / a[n][n]
for k in range(n-1, 0, -1):
akjxj = 0.0
for j in range(k+1, n+1):
akjxj += a[k][j] * x[j]
x[k] = (b[k] - akjxj) / a[k][k]
def main():
# 配列は0番目ダミー
f = [0.0] * 19 #節点力ベクトル
d = [0.0] * 19 #節点変位ベクトル
x = [0.0] * 19 #境界条件（拘束節点に0, その他に1が入る)
ss = [[0.0]*19 for _ in range(19)] #全体剛性行列
xv = [0.0] * 10 #節点ごとの変位境界 y方向
xw = [0.0] * 10 #節点ごとの変位境界 z方向
fy = [0.0] * 10 #節点ごとの荷重条件 y方向
fz = [0.0] * 10 #節点ごとの荷重条件 z方向
hidari = [0] * 10 #各要素の左節点番号
migi = [0] * 10 #各要素の右節点番号
s = [[[0.0]*5 for _ in range(5)] for _ in range(10)] #要素剛性行列
th = [0.0] * 10 #各要素の回転角
t = [[0.0]*5 for _ in range(5)] #座標変換行列（1要素ごとに計算）
s44 = [[0.0]*5 for _ in range(5)] #要素剛性行列（1要素ごとに移し替える入れ物）
ts = [[0.0]*5 for _ in range(5)] #TK（1要素ごとに計算）
tst = [[0.0]*5 for _ in range(5)] #TKT（1要素ごとに計算）
sn = [0.0] * 5 #節点力の出力のために追加
pi = math.pi
nyou = 2 #要素数
nset = 3 #節点数
# 初期化
for n in range(1, nyou+1):
for i in range(1, 5):
for j in range(1, 5):
s[n][i][j] = 0.0
for i in range(1, 19):
d[i] = 0.0
# f[i] = 0.0
# x[i] = 0.0
for j in range(1, 19):
ss[i][j] = 0.0
for i in range(1, 10):
xv[i] = 1.0
xw[i] = 1.0
fy[i] = 0.0
fz[i] = 0.0
E = 1.0 #ヤング率 N/mm^2(MPa)
ell = 1.0 #1要素の長さ mm
A = 1.0 #面積 mm^2
sk1 = E*A/ell #要素1のばね定数
# 要素剛性マトリクス
s[1][2][2] = sk1; s[1][2][4] = -sk1
s[1][4][2] = -sk1; s[1][4][4] = sk1
sk2 = sk1 / math.sqrt(2.0) #要素2のばね定数
s[2][2][2] = sk2; s[2][2][4] = -sk2
s[2][4][2] = -sk2; s[2][4][4] = sk2
#各要素の左節点番号、右節点番号
hidari[1] = 1; migi[1] = 2
hidari[2] = 2; migi[2] = 3
#各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th[1] = 0.0
th[2] = pi + pi/4.0
# 境界条件
xv[1] = 0.0; xw[1] = 0.0
xv[3] = 0.0; xw[3] = 0.0
for i in range(1, 10):
x[2*i-1] = xv[i]
x[2*i] = xw[i]
# 荷重条件(y方向:fy[節点番号], z方向:fz[節点番号])
fy[2] = 1.0 #[N]
for i in range(1, 10):
f[2*i-1] = fy[i]
f[2*i] = fz[i]
# 要素ごとのTKTの計算
for n in range(1, nyou+1):
zahyou(t, th[n])
for i in range(1, 5):
for j in range(1, 5):
s44[i][j] = s[n][i][j]
mxtmx(t, s44, ts)
mxmx(ts, t, tst)
i1 = 2*hidari[n] - 1
i2 = 2*hidari[n]
m1 = 2*migi[n] - 1
m2 = 2*migi[n]
ss[i1][i1] += tst[1][1]
ss[i1][i2] += tst[1][2]
ss[i1][m1] += tst[1][3]
ss[i1][m2] += tst[1][4]
ss[i2][i1] += tst[2][1]
ss[i2][i2] += tst[2][2]
ss[i2][m1] += tst[2][3]
ss[i2][m2] += tst[2][4]
ss[m1][i1] += tst[3][1]
ss[m1][i2] += tst[3][2]
ss[m1][m1] += tst[3][3]
ss[m1][m2] += tst[3][4]
ss[m2][i1] += tst[4][1]
ss[m2][i2] += tst[4][2]
ss[m2][m1] += tst[4][3]
ss[m2][m2] += tst[4][4]
# 境界条件を入れる
for i in range(1, nset*2+1):
for j in range(1, nset*2+1):
ss[i][j] = x[i]*ss[i][j]
ss[j][i] = x[i]*ss[j][i]
for i in range(1, nset*2+1):
if x[i] < 1e-3:
ss[i][i] = 1.0
gausu(nset*2, ss, d, f)
for n in range(1, nset+1):
print(f"節点番号：{n} v={d[n*2-1]:.10f}, w={d[n*2]:.10f}")
# 要素ごとのTKTの計算
for n in range(1, nyou+1):
zahyou(t, th[n])
for i in range(1, 5):
for j in range(1, 5):
s44[i][j] = s[n][i][j]
mxtmx(t, s44, ts)
mxmx(ts, t, tst)
for i in range(1, 5):
sn[i] = 0.0
for j in range(1, 3):
sn[i] += tst[i][j] * d[(hidari[n]-1)*2 + j]
for j in range(3, 5):
sn[i] += tst[i][j] * d[(migi[n]-1)*2 + j-2]
print(f"部材番号： {n} S({hidari[n]})={sn[1]:15.10f}, N({hidari[n]})={sn[2]:15.10f}")
print(f" S({migi[n]})={sn[3]:15.10f}, N({migi[n]})={sn[4]:15.10f}")
if __name__ == "__main__":
main()