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.asin(1.0) * 2.0 nyou = 5 #要素数 nset = 4 #節点数 # 初期化 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 = 1000.0 #1要素の長さ[mm] ea = 10.0e3 * 10.0**2 #EA, ヤング率[N/mm^2]*面積[mm^2] ei = 10.0e3 * 10.0**4 / 12.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から要素4までの剛性マトリクス for k in range(2, 5): for i in range(1, 7): for j in range(1, 7): s[k][i][j] = s[1][i][j] #要素5の剛性マトリクス el = el * math.sqrt(2.0) s[5][2][2] = ea/el s[5][2][5] = -ea/el s[5][5][2] = -ea/el s[5][5][5] = ea/el s[5][1][1] = 12.0*ei/el**3 s[5][1][3] = -6.0*ei/el**2 s[5][1][4] = -12.0*ei/el**3 s[5][1][6] = -6.0*ei/el**2 s[5][3][1] = -6.0*ei/el**2 s[5][3][3] = 4.0*ei/el s[5][3][4] = 6.0*ei/el**2 s[5][3][6] = 2.0*ei/el s[5][4][1] = -12.0*ei/el**3 s[5][4][3] = 6.0*ei/el**2 s[5][4][4] = 12.0*ei/el**3 s[5][4][6] = 6.0*ei/el**2 s[5][6][1] = -6.0*ei/el**2 s[5][6][3] = 2.0*ei/el s[5][6][4] = 6.0*ei/el**2 s[5][6][6] = 4.0*ei/el #各要素の左節点番号、右節点番号 hidari[1]=1; migi[1]=2 hidari[2]=2; migi[2]=3 hidari[3]=3; migi[3]=4 hidari[4]=4; migi[4]=1 hidari[5]=1; migi[5]=3 #各要素の回転角（上の左右節点番号がz軸に横たわる状態からの） th[1]=0.0 th[2]=pi/2.0 th[3]=pi th[4]=pi+pi/2.0 th[5]=pi/4.0 # 境界条件 xv[1]=0.0 xw[1]=0.0 xv[2]=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]=500.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}")