マトリクス構造解析オンライン授業用テキスト
第14回オンライン授業
YouTube授業動画リスト
小さい字は補足説明なので、読み飛ばしてもいいです。
マトリクス構造解析Iの第8回で、
片持ち折れ梁を単位荷重法で解いたが、今回は、
これをプログラムで解いてみたい。
単位荷重法で求めた載荷点の変位は、
$v_{3}
=-\frac{P\ell^{3}}{2EI}$
$w_{3}=\frac{P\ell}{EA}+\frac{4P\ell^{3}}{3EI}$
で、断面力は、
$S_{①}=0$
$N_{①}=P$
$M_{①}=P\ell$
$S_{②}=-P$
$N_{②}=0$
$M_{②}=P(\ell-z_{②}^{\ell})$ 部材端での値は、
$M_{②}(0)=P\ell, \;M_{②}(\ell)=0$
部材①、②を切り離して、各節点での断面力の値を、
正の値で描ける向きに矢印を向けて描くと図のようになる。
つまり、各節点の節点力として$y, z$方向、$x$軸右ねじ方向プラスで書くと
以下のようになる。
全体系での節点力なので、
$y$方向であれば軸力だろうがせん断力だろうが$S$となり、
$z$方向であれば軸力だろうがせん断力だろうが$N$となることに注意すると、
以下のようになる。
$S_{1}^{①}=0$
$N_{1}^{①}=-P$
$M_{1}^{①}=-P\ell$
$S_{2}^{①}=0$
$N_{2}^{①}=P$
$M_{2}^{①}=P\ell$
$S_{2}^{②}=0$
$N_{2}^{②}=-P$
$M_{2}^{②}=-P\ell$
$S_{3}^{②}=0$
$N_{3}^{②}=P$
$M_{3}^{②}=0$
となる。
では、 梁の剛性方程式を用いた 骨組のプログラムを使って、上の折れ梁を解いてみよう。 まず、 hone2.pyを ダウンロードするか、 このページの末尾のソースをテキストエディタに貼り付けて保存する。 使い方は、torasu2.pyとほぼ同じだが、一応、説明する。
今回は、要素数は2要素、節点数は3節点。
nyou = 2 #要素数 nset = 3 #節点数
まず、以下のように書いてあるところを見つける。
# 要素1の剛性マトリクス el = 1.0 #1要素の長さ[mm] ea = 1.0 #EA, ヤング率[N/mm^2]*面積[mm^2] ei = 1.0 #EI, ヤング率[N/mm^2]*断面2次モーメント[mm^4] s[1][2][2] = ea/el s[1][2][5] = -ea/el s[1][5][2] = -ea/el s[1][5][5] = ea/el s[1][1][1] = 12.0*ei/el**3 s[1][1][3] = -6.0*ei/el**2 s[1][1][4] = -12.0*ei/el**3 s[1][1][6] = -6.0*ei/el**2 s[1][3][1] = -6.0*ei/el**2 s[1][3][3] = 4.0*ei/el s[1][3][4] = 6.0*ei/el**2 s[1][3][6] = 2.0*ei/el s[1][4][1] = -12.0*ei/el**3 s[1][4][3] = 6.0*ei/el**2 s[1][4][4] = 12.0*ei/el**3 s[1][4][6] = 6.0*ei/el**2 s[1][6][1] = -6.0*ei/el**2 s[1][6][3] = 2.0*ei/el s[1][6][4] = 6.0*ei/el**2 s[1][6][6] = 4.0*ei/el
torasu2.pyと同様に、 要素nの要素剛性マトリクス$k_{ij}$が、s(n,i,j)として 与えられている。 これにマトリクス構造解析Iの第7回でやった軸力ありの梁の剛性方程式を入れてやればよい。 上記の要素①については、既に、$6\times 6$のマトリクスの成分が、 ea($EA$), el($\ell$), ei($EI$)などの変数で与えられているから、 あとは、ea, el, eiを適切に与えてやればよい。 マトリクスは初期化しているので、0の成分は与えなくてもよい。 今回は、第4回でやったのと同様に、 計算を確かめるために、$EA=1, EI=1, \ell=1$を与えて計算する。 つまり、el, ea, eiのところをすべて1.0に書き換えればよい。
el = 1.0 #1要素の長さ[mm] ea = 1.0 #EA, ヤング率[N/mm^2]*面積[mm^2] ei = 1.0 #EI, ヤング率[N/mm^2]*断面2次モーメント[mm^4]
次にその下の要素②以降のマトリクスを定義する。
#要素2の剛性マトリクス
for i in range(1, 7):
for j in range(1, 7):
s[2][i][j] = s[1][i][j]
上の例は、要素②の剛性マトリクスs[2][i][j] (つまり$k_{ij}^{②}$) に要素①の剛性マトリクスs[1][i][j] (つまり$k_{ij}^{①}$)を代入している。 添字i,jは、1から6まで変化させればいいが、 for文の初期値は1からでも回数は0から数えるので7回で range(1,7)となる。
今回の要素は、初期状態で$z$軸に左から要素①、②、節点1, 2, 3の
順で並べられるから、
この状態での各要素の左節点番号(hidari)と右節点番号(migi)を与える。
要素③以降は、削除するか ! を入れてコメントにする。
#各要素の左節点番号、右節点番号 hidari[1]=1; migi[1]=2 hidari[2]=2; migi[2]=3
次に各要素の回転角を、初期状態での左節点を中心に反時計回りの回転角で与える。 $yz$平面の回転角は、$x$軸右ねじ回転なので、反時計回りになる。 要素①は初期状態から動かないのでth(1)=0. 要素②は、$\pi+\frac{\pi}{2}$だけ 回転させればよい。 要素③以降はないので、削除するか、! を入れてコメントにする。
#各要素の回転角（上の左右節点番号がz軸に横たわる状態からの） th[1]=0.0 th[2]=pi+pi/2.0
次は境界条件である。 回転角の下の以下のように書いてあるところを見つける。
# 境界条件 xv[1]=0.0 xw[1]=0.0 xt[1]=0.0
[]内は拘束する節点番号、xvは$v$, xwは$w$, xtは$\theta$を意味する。 今回は、節点1で、$v_{1}=w_{1}=\theta_{1}=0$だから、 上のように与える。
境界条件の下の 以下のように書いてあるとこを見つける。
# 荷重条件(y方向:fy[節点番号], z方向:fz[節点番号], モーメント外力:cx[節点番号]) fz[3]=1.0 #[N]
[]内は、荷重を与える節点番号。 fyは$y$方向荷重、fzは$z$方向荷重、cxは$x$軸右ねじ回りのモーメント荷重を与える。 今回は、節点3の$z$方向に$P=1$を与えたいので、上のようにする。
プログラムの修正が終わったら、 修正したプログラムは、例えば"hone2.py"みたいに名前を変えて保存する。 プログラムを実行すると、 以下のような結果が表示される。
節点番号：1 v=0.0, w=0.0, th=0.0
節点番号：2 v=-0.5000000000000001, w=0.9999999999999998, th=1.0000000000000002
節点番号：3 v=-0.49999999999999994, w=2.333333333333333, th=1.5
曲げモーメントの単位に注意!
部材番号： 1 S(1)= 0.0000000000, N(1)= -1.0000000000, M(1)= -1.0000000000
S(2)= 0.0000000000, N(2)= 1.0000000000, M(2)= 1.0000000000
部材番号： 2 S(2)= -0.0000000000, N(2)= -1.0000000000, M(2)= -1.0000000000
S(3)= 0.0000000000, N(3)= 1.0000000000, M(3)= 0.0000000000
単位荷重法では、
$v_{3}
=-\frac{P\ell^{3}}{2EI}$
$w_{3}=\frac{P\ell}{EA}+\frac{4P\ell^{3}}{3EI}$
だから、$P=1, \ell=1, EI=1, EA=1$を代入すると、
$v_{3}=-\frac{1}{2}=-0.5$
$w_{3}=1+\frac{4}{3}=2.333333....$
であるから、節点番号3のところを見ると、ほぼ合っている。
節点力は、$P=1, \ell=1$とすると、
$
\begin{array}{lll}
S_{1}^{①}=0, &
N_{1}^{①}=-P=-1, &
M_{1}^{①}=-P\ell=-1 \\
S_{2}^{①}=0, &
N_{2}^{①}=P=1, &
M_{2}^{①}=P\ell=1 \\
S_{2}^{②}=0, &
N_{2}^{②}=-P=-1, &
M_{2}^{②}=-P\ell=-1 \\
S_{3}^{②}=0, &
N_{3}^{②}=P=1, &
M_{3}^{②}=0 \\
\end{array}
$
だから、これも合ってそうだ。
以下の片持ち折れ梁（自分の学籍番号のもの）の 載荷点の節点変位($v, w$)を単位荷重法で求めよ。 部材の長さ、伸び剛性、曲げ剛性は2部材ともそれぞれ$\ell$, $EA$, $EI$とする。 次に、$P=E=A=\ell=1.$として、 それをhone2.pyで解いて、 手計算の答えとどれくらい合うか考察せよ。 WebClassの「プログラム課題6」から、 単位荷重法の手計算を撮影した画像と プログラム本体（例えば hone2.py とか）と、 出力ファイル（例えば kekka6.txt とか）を それぞれアップロードせよ。 どうしてもエラーが取れないなどの場合は、 kekka6.txtの代わりに、エラー画面等 をテキストエディタに貼り付けて保存したものをアップロードする。 また、最後のコメント欄に、手計算とプログラムの計算とを比較した考察や、 うまく実行できなかった場合は、 どういう症状で実行できないのかの説明を書き込むこと。
学籍番号の末尾が1か8の人。
学籍番号の末尾が2か9の人。
学籍番号の末尾が3か0の人。
学籍番号の末尾が4か6の人。
学籍番号の末尾が5か7の人。
import math
# 1始まり配列を作るユーティリティ
def zeros1d(n):
return [0.0] * (n + 1) # 0番目はダミー
def zeros2d(n, m):
return [[0.0] * (m + 1) for _ in range(n + 1)]
def zeros3d(n, m, l):
return [[[0.0] * (l + 1) for _ in range(m + 1)] for _ in range(n + 1)]
# ---- サブルーチン ----
def zahyou(t, theta):
for i in range(1, 4):
for j in range(1, 4):
t[i+3][j] = 0.0
t[i][j+3] = 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] = 1.0
t[4][4] = t[1][1]
t[4][5] = t[1][2]
t[5][4] = t[2][1]
t[5][5] = t[2][2]
t[6][6] = 1.0
def mxmx(a, b, c):
for i in range(1, 7):
for j in range(1, 7):
c[i][j] = 0.0
for k in range(1, 7):
c[i][j] += a[i][k] * b[k][j]
def mxtmx(a, b, c):
for i in range(1, 7):
for j in range(1, 7):
c[i][j] = 0.0
for k in range(1, 7):
c[i][j] += a[k][i] * b[k][j]
def gausu(n, a, x, b):
# 前進消去
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, 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]
# ---- メインプログラム ----
# 配列定義
f = zeros1d(27) #節点力ベクトル
d = zeros1d(27) #節点変位ベクトル
x = zeros1d(27) #境界条件（拘束節点に0, その他に1が入る)
ss = zeros2d(27, 27) #全体剛性行列
xv = zeros1d(9) #節点ごとの変位境界 y方向
xw = zeros1d(9) #節点ごとの変位境界 z方向
xt = zeros1d(9) #節点ごとの変位境界 回転
fy = zeros1d(9) #節点ごとの荷重条件 z方向
fz = zeros1d(9) #節点ごとの荷重条件 y方向
cx = zeros1d(9) #節点ごとの荷重条件 モーメント外力
hidari = zeros1d(9) #各要素の左節点番号
migi = zeros1d(9) #各要素の右節点番号
s = zeros3d(9, 6, 6) #要素剛性行列
th = zeros1d(9) #各要素の回転角
t = zeros2d(6, 6) #座標変換行列（1要素ごとに計算）
s66 = zeros2d(6, 6) #要素剛性行列（1要素ごとに移し替える入れ物）
ts = zeros2d(6, 6) #TK（1要素ごとに計算）
tst = zeros2d(6, 6) #TKT（1要素ごとに計算）
snm = zeros1d(6) #節点力の出力のために追加
pi = math.pi
nyou = 2 #要素数
nset = 3 #節点数
# 初期化
for n in range(1, nyou+1):
for i in range(1, 7):
for j in range(1, 7):
s[n][i][j] = 0.0
for i in range(1, 28):
d[i] = 0.0
for j in range(1, 28):
ss[i][j] = 0.0
for i in range(1, 10):
xv[i] = 1.0
xw[i] = 1.0
xt[i] = 1.0
fy[i] = 0.0
fz[i] = 0.0
cx[i] = 0.0
# 要素1の剛性マトリクス
el = 1.0 #1要素の長さ[mm]
ea = 1.0 #EA, ヤング率[N/mm^2]*面積[mm^2]
ei = 1.0 #EI, ヤング率[N/mm^2]*断面2次モーメント[mm^4]
s[1][2][2] = ea/el
s[1][2][5] = -ea/el
s[1][5][2] = -ea/el
s[1][5][5] = ea/el
s[1][1][1] = 12.0*ei/el**3
s[1][1][3] = -6.0*ei/el**2
s[1][1][4] = -12.0*ei/el**3
s[1][1][6] = -6.0*ei/el**2
s[1][3][1] = -6.0*ei/el**2
s[1][3][3] = 4.0*ei/el
s[1][3][4] = 6.0*ei/el**2
s[1][3][6] = 2.0*ei/el
s[1][4][1] = -12.0*ei/el**3
s[1][4][3] = 6.0*ei/el**2
s[1][4][4] = 12.0*ei/el**3
s[1][4][6] = 6.0*ei/el**2
s[1][6][1] = -6.0*ei/el**2
s[1][6][3] = 2.0*ei/el
s[1][6][4] = 6.0*ei/el**2
s[1][6][6] = 4.0*ei/el
#要素2の剛性マトリクス
for i in range(1, 7):
for j in range(1, 7):
s[2][i][j] = s[1][i][j]
#各要素の左節点番号、右節点番号
hidari[1]=1; migi[1]=2
hidari[2]=2; migi[2]=3
#各要素の回転角（上の左右節点番号がz軸に横たわる状態からの）
th[1]=0.0
th[2]=pi+pi/2.0
# 境界条件
xv[1]=0.0
xw[1]=0.0
xt[1]=0.0
for i in range(1, 10):
x[3*i-2] = xv[i]
x[3*i-1] = xw[i]
x[3*i] = xt[i]
# 荷重条件(y方向:fy[節点番号], z方向:fz[節点番号], モーメント外力:cx[節点番号])
fz[3]=1.0 #[N]
for i in range(1, 10):
f[3*i-2] = fy[i]
f[3*i-1] = fz[i]
f[3*i] = cx[i]
# 要素ごとのTKTの計算
for n in range(1, nyou+1):
zahyou(t, th[n])
for i in range(1, 7):
for j in range(1, 7):
s66[i][j] = s[n][i][j]
mxtmx(t, s66, ts)
mxmx(ts, t, tst)
i1 = 3*hidari[n]-2
i2 = 3*hidari[n]-1
i3 = 3*hidari[n]
m1 = 3*migi[n]-2
m2 = 3*migi[n]-1
m3 = 3*migi[n]
ss[i1][i1] += tst[1][1]
ss[i1][i2] += tst[1][2]
ss[i1][i3] += tst[1][3]
ss[i1][m1] += tst[1][4]
ss[i1][m2] += tst[1][5]
ss[i1][m3] += tst[1][6]
ss[i2][i1] += tst[2][1]
ss[i2][i2] += tst[2][2]
ss[i2][i3] += tst[2][3]
ss[i2][m1] += tst[2][4]
ss[i2][m2] += tst[2][5]
ss[i2][m3] += tst[2][6]
ss[i3][i1] += tst[3][1]
ss[i3][i2] += tst[3][2]
ss[i3][i3] += tst[3][3]
ss[i3][m1] += tst[3][4]
ss[i3][m2] += tst[3][5]
ss[i3][m3] += tst[3][6]
ss[m1][i1] += tst[4][1]
ss[m1][i2] += tst[4][2]
ss[m1][i3] += tst[4][3]
ss[m1][m1] += tst[4][4]
ss[m1][m2] += tst[4][5]
ss[m1][m3] += tst[4][6]
ss[m2][i1] += tst[5][1]
ss[m2][i2] += tst[5][2]
ss[m2][i3] += tst[5][3]
ss[m2][m1] += tst[5][4]
ss[m2][m2] += tst[5][5]
ss[m2][m3] += tst[5][6]
ss[m3][i1] += tst[6][1]
ss[m3][i2] += tst[6][2]
ss[m3][i3] += tst[6][3]
ss[m3][m1] += tst[6][4]
ss[m3][m2] += tst[6][5]
ss[m3][m3] += tst[6][6]
# 境界条件を入れる
for i in range(1, nset*3+1):
for j in range(1, nset*3+1):
ss[i][j] = x[i] * ss[i][j]
ss[j][i] = x[i] * ss[j][i]
for i in range(1, nset*3+1):
if x[i] < 1e-3:
ss[i][i] = 1.0
# 線形方程式を解く
gausu(nset*3, ss, d, f)
# 結果出力
for n in range(1, nset+1):
print(f"節点番号：{n} v={d[n*3-2]}, w={d[n*3-1]}, th={d[n*3]}")
print(" 曲げモーメントの単位に注意!")
for n in range(1, nyou+1):
zahyou(t, th[n])
for i in range(1, 7):
for j in range(1, 7):
s66[i][j] = s[n][i][j]
mxtmx(t, s66, ts)
mxmx(ts, t, tst)
for i in range(1, 7):
snm[i] = 0.0
for j in range(1, 4):
snm[i] += tst[i][j] * d[(hidari[n]-1)*3 + j]
for j in range(4, 7):
snm[i] += tst[i][j] * d[(migi[n]-1)*3 + j-3]
print(f"部材番号： {n} S({hidari[n]})={snm[1]:15.10f}, N({hidari[n]})={snm[2]:15.10f}, M({hidari[n]})={snm[3]:15.10f}")
print(f" S({migi[n]})={snm[4]:15.10f}, N({migi[n]})={snm[5]:15.10f}, M({migi[n]})={snm[6]:15.10f}")