マトリクス変位法ではりの計算式を解くやり方を紹介しています。はりの中間の応力と変位を弾性曲線式によって解いています。
本稿の弾性曲線式は
V''''(x) = p(x), x | 0 => L (technical theory of beam, Euler–Bernoulli beam theory)
V''''(x) = p(x) - κ * EI / GA * p''(x), x | 0 => L (Timoshenko beam theory)
と改めます。上は普通の梁で、下はせん断変形を考慮した梁の式です。
はりの計算式は
V''''(x) = p(x) = -w (1)
V''' (x) = Q(x) = C1 + ∫ p(x) dx (2)
V'' (x) = M(x) = C2 + ∫ Q(x) dx (3)
r(x) = -M(x) (4)
G(x) = κ * EI / GA * Q(x) (4-1)
V' (x) = T(x) = C3 + ∫ r(x) dx (5)
V (x) = V(x) = C4 + ∫ T(x) + G(x) dx (6)
ただし
κ せん断変形の形状係数
EI 曲げ剛性(一定)
GA せん断剛性(一定)
x 位置
L スパン長
p(x) 分布荷重
w 等分布荷重
Q(x) せん断力
M(x) 曲げモーメント
r(x) 曲率(δ''=φ=r(x)/EI)
G(x) せん断角(γ=G(x)/EI、鉛直方向のせん断歪は全断面で一定とし、材軸方向は０とする)
T(x) たわみ角(θ=T(x)/EI)
V(x) たわみ(δ=V(x)/EI)
C1 積分定数(Q(x=0))
C2 積分定数(M(x=0))
C3 積分定数(T(x=0))
C4 積分定数(V(x=0))
とします。G(x) = 0 とすれば普通の梁です。
台形分布荷重で、はりの計算式をあらかじめ積分すると
pl = wf * d (5-1)
Q(x) = C1 +
-w1 / 2 * (x - a) - w2 * (x - a - d) - w3 / 2 * (x - a - 2 * b + d) +
-p1 - p2 - p3 (5-2)
M(x) = C2 + C1 * x +
-w1 / 6 * (x - a)^2 - w2 / 2 * (x - a - d)^2 - w3 / 6 * (x - a - 2 * b + d)^2 +
-p1 * (x - a - 2 / 3 * d) - p2 * (x - a - d - (2 * b - d) / 2) - p3 * (x - a - 2 * b + d - 2 / 3 * d) (5-3)
G(x) =(C1 +
-w1 / 2 * (x - a) - w2 * (x - a - d) - w3 / 2 * (x - a - 2 * b + d) +
-p1 - p2 - p3) * __g__ * L^2 (5-4)
T(x) = C3 - C2 * x - C1 / 2.0 * x^2 +
w1 / 24 * (x - a)^3 + w2 / 6 * (x - a - d)^3 + w3 / 24 * (x - a - 2 * b + d)^3 +
p1 * d^2 / 12 + p2 * (2 * b - d)^2 / 6 + p3 * d^2 / 12 +
p1 * d / 3 * (x - a - d) + p2 * (2 * b - d) / 2 * (x - a - 2 * b) + p3 * d / 3 * (x - a - 2 * b) +
p1 / 2 * (x - a - d)^2 + p2 / 2 * (x - a - 2 * b)^2 + p3 / 2 * (x - a - 2 * b)^2 (5-5)
V(x) = C4 + C3 * x - C2 / 2 * x^2 - C1 / 6 * x^3 +
w1 / 120 * (x - a)^4 + w2 / 24 * (x - a - d)^4 + w3 / 120 * (x - a - 2 * b + d)^4 +
p1 * d^2 / 12 * (x - a - d * 4 / 5) +
p2 * (2 * b - d)^2 / 6 * (x - a - d - (2 * b - d) * 3 / 4) +
p3 * d^2 / 12 * (x - a - 2 * b + d - d * 4 / 5) +
p1 * d / 3 * (x - a - d)^2 / 2 +
p2 * (2 * b - d) / 2 * (x - a - 2.0 * b)^2 / 2 +
p3 * d / 3 * (x - a - 2 * b)^2 / 2 +
p1 / 6 * (x - a - d)^3 + p2 / 6 * (x - a - 2 * b)^3 + p3 / 6 * (x - a - 2 * b)^3 + __g__ * (m - ma) * L^2 (5-6)
ただし
__g__ = κ * EI / GA / L^2 ティモシェンコ数
b = (l - a - c) / 2
w1 = pl * ( (x - a) / d) * (x < a ? 0 : (x < a + d ? 1 : 0))
w2 = pl * (x < a + d ? 0 : (x < a + 2 * b ? 1 : 0))
w3 =-pl * ( (x - a - 2 * b + d) / d) * (x < a + 2 * b - d ? 0 : (x < a + 2 * b ? 1 : 0))
p1 = pl * d / 2 * (x < a + d ? 0 : 1)
p2 = pl * (2 * b - d) * (x < a + 2 * b ? 0 : 1)
p3 =-pl * d / 2 * (x < a + 2 * b ? 0 : 1)
となります。__g__ = 0 とすれば普通の梁です。
はり端の Q(x=0)、M(x=0)、T(x=0)、V(x=0) が積分定数となるので、マトリクス変位法で算定して、積分定数を代入すれば計算式は解けます。
Mery で以下のスクリプトを ruby で書いて ファイル名 "000-16.bat" で保存して実行してください。
１端ピン他端固定、台形分布荷重としています。
: マトリクス変位法(台形分布荷重)
@echo off
ruby -x %~f0 4 5 1 1 1 -10 0.15 0.3 6.5 6.5/15 1.2
pause
goto:eof
#!ruby -rmatrix
# 境界条件
def boundary_condition(tp, l, ca, cb, qa, qb, ei = nil, __g__ = 0, m1 = 0, m2 = 0, ka = 0, kb = 0, k1 = 1e20, k2 = 1e20)
=begin
tp 支持条件
0 単純支持
1 Ａ端自由・Ｂ端固定
2 Ａ端固定・Ｂ端自由
3 両端固定
4 Ａ端ピン・Ｂ端固定
5 Ａ端固定・Ｂ端ピン
6 両端回転・鉛直バネ支持
l スパン m
ca Ａ端の固定端モーメント kNm
cb Ｂ端の固定端モーメント kNm
qa Ａ端の単純支持のせん断力 kN
qb Ｂ端の単純支持のせん断力 kN
ei 曲げ剛性 kN*m2
__g__ ティモシェンコ数
m1 Ａ端の回転外力 kNm
m2 Ｂ端の回転外力 kNm
ka Ａ端支持の回転バネ kNm/rad
kb Ｂ端支持の回転バネ kNm/rad
k1 Ａ端支持の鉛直バネ kN/m
k2 Ｂ端支持の鉛直バネ kN/m
( [k] + [d] ) * {y} + {c} = {f}
[k] : 部材剛性マトリクス
[d] : 支点バネマトリクス
{y} : 節点の変位ベクトル → [ va, ta, vb, tb ] = ( [k] + [d] )^^-1 * ( {f} - {c} )
{c} : 固定材端力ベクトル
{f} : 節点の荷重ベクトル( 節点外力 )
{r} = [k] * {y} + {c}
{r} : 材端の内力ベクトル → [-ra, ma,-rb, mb ] = [k] * {y} + {c}
ra Ａ端のせん断力 kN
ma Ａ端の曲げモーメント kNm
ta Ａ端のたわみ角 rad / ei
va Ａ端のたわみ m / ei
rb Ｂ端のせん断力 kN
mb Ｂ端の曲げモーメント kNm
tb Ｂ端のたわみ角 rad / ei
vb Ｂ端のたわみ m / ei
=end
# A) 剛性マトリクス [k]
a1 = 12.0 / l ** 3
b1 = 6.0 / l ** 2
c1 = b1
d1 = 2.0 / l
f1 = 4.0 / l
h1 = f1
if 0 < __g__ && ei
# せん断変形
p3 = 12.0 * __g__
a1 *= 1 / (1 + p3)
b1 *= 1 / (1 + p3)
c1 *= 1 / (1 + p3)
d1 *= 1 / (1 + p3) * (1 - p3 / 2)
f1 *= 1 / (1 + p3) * (1 + p3 / 4)
h1 *= 1 / (1 + p3) * (1 + p3 / 4)
end
k = Matrix[
[ a1, b1,-a1, c1],
[ b1, f1,-b1, d1],
[-a1,-b1, a1,-c1],
[ c1, d1,-c1, h1] ].map { |u| u.to_f }
# B) 固定材端力ベクトル {c}
c = Matrix[
[-da = qa - (ca + cb) / l],
[ ca],
[-db = qb + (ca + cb) / l],
[ cb] ].map { |u| u.to_f }
# C) 節点外力 {f} = {f} - {c}
f = Matrix[
[ p1 = 0],
[ m1 = m1],
[ p2 = 0],
[ m2 = m2] ].map { |u| u.to_f } - c
# D) 支点バネ [d]
case tp
when 0; k1 = k2 = 1e20; ka = kb = 0 # 3-1) 単純支持
when 1; k1 = ka = 0; k2 = kb = 1e20 # 3-2) Ａ端自由・Ｂ端固定
when 2; k2 = kb = 0; k1 = ka = 1e20 # 3-3) Ａ端固定・Ｂ端自由
when 3; k1 = k2 = ka = kb = 1e20 # 3-4) 両端固定
when 4; k1 = k2 = kb = 1e20; ka = 0 # 3-5) Ａ端ピン・Ｂ端固定
when 5; k1 = k2 = ka = 1e20; kb = 0 # 3-6) Ａ端固定・Ｂ端ピン
when 6 # 3-7) 両端回転・鉛直バネ支持
else
k1 = k2 = 1e20; ka = kb = 0
end
d = Matrix[
[ k1, 0, 0, 0],
[ 0, ka, 0, 0],
[ 0, 0, k2, 0],
[ 0, 0, 0, kb] ].map { |u| u.to_f }
# E) 節点の変位ベクトル {y} = ( [k] + [d] )^^-1 * {f}
y = (k + d / ei).inv * f
# F) 材端の内力ベクトル {r} = [k] * {y} + {c}
r = k * y + c
va, ta, vb, tb = y.to_a.flatten
ra, ma, rb, mb = r.to_a.flatten
return -ra,-rb, ma, mb, ta, tb, va, vb
end
def cq5(wf, l, a, c, d, __g__ = 0)
=begin
w = wf * d * (2 * b - d)
w = wf * d * (l - a - c - d)
|d|
A _____wf B ___
_____／_____＼_____ _|_d
^ ^
|-a-|---2*b---|-c-|
0 L ---> x
|--------L--------|
l スパン m
wf 荷重 kN/m2
a 荷重の開始位置 m
c 荷重の終了位置 m
d 荷重の奥行 m
__g__ ティモシェンコ数
ca Ａ端の固定端モーメント kNm
cb Ｂ端の固定端モーメント kNm
qa Ａ端の単純支持のせん断力 kN
qb Ｂ端の単純支持のせん断力 kN
=end
b = (l - a - c) / 2.0
w = wf * d * (2 * b - d)
ca =-w / (6 * l ** 2) * ( (b ** 2 + (d - b) ** 2 + 6 * (c + b) ** 2) * (a - b - 2 * c) + 12 * (c + b) ** 3)
cb = w / (6 * l ** 2) * ( (b ** 2 + (d - b) ** 2 + 6 * (a + b) ** 2) * (c - b - 2 * a) + 12 * (a + b) ** 3)
qa = w / l * (c + b)
qb = w / l * (a + b)
if __g__ != 0
cs = w / (6 * l ** 2) * (6 * (a + b) * (c + b) - b ** 2 - (d - b) ** 2) * 6 * __g__ * l
ca = (ca - cs) / (12 * __g__ + 1)
cb = (cb + cs) / (12 * __g__ + 1)
end
return ca, cb, qa, qb
end
def mq5(x, l, a, c, d, wf, ra, ma, ta, va, ei = nil, __g__ = 0, ka = 0, kb = 0, k1 = 1e20, k2 = 1e20)
=begin
x 位置 m
l スパン m
wf 荷重 kN/m2
a 荷重の開始位置 m
c 荷重の終了位置 m
d 荷重の奥行 m
ra Ａ端のせん断力 kN
ma Ａ端の曲げモーメント kNm
ta Ａ端のたわみ角 rad
va Ａ端のたわみ m
ei 曲げ剛性 kN*m2
__g__ ティモシェンコ数
|d|
A _____wf B ___
_____／_____＼_____ _|_d
^ ^
|-a-|---2*b---|-c-|
0 L ---> x
|--------L--------|
_____
____ ／_____＼ ____
1____ ／|______ ____
________
2____ __| |____
3____ _________ ____
＼|
q せん断力 kN
m 曲げモーメント kNm
t / ei たわみ角 rad
v / ei たわみ m
=end
b = (l - a - c) / 2.0
pl = wf * d
w1 = pl * ( (x - a) / d) * (x < a ? 0 : (x < a + d ? 1 : 0))
w2 = pl * (x < a + d ? 0 : (x < a + 2.0 * b ? 1 : 0))
w3 =-pl * ( (x - a - 2.0 * b + d) / d) * (x < a + 2.0 * b - d ? 0 : (x < a + 2.0 * b ? 1 : 0))
p1 = pl * d / 2.0 * (x < a + d ? 0 : 1)
p2 = pl * (2.0 * b - d) * (x < a + 2.0 * b ? 0 : 1)
p3 =-pl * d / 2.0 * (x < a + 2.0 * b ? 0 : 1)
q = ra +
-w1 / 2.0 * (x - a) - w2 * (x - a - d) - w3 / 2.0 * (x - a - 2.0 * b + d) +
-p1 - p2 - p3
m = ma + ra * x +
-w1 / 6.0 * (x - a) ** 2 - w2 / 2.0 * (x - a - d) ** 2 - w3 / 6.0 * (x - a - 2.0 * b + d) ** 2 +
-p1 * (x - a - 2.0 / 3.0 * d) - p2 * (x - a - d - 1.0 / 2.0 * (2.0 * b - d)) - p3 * (x - a - 2.0 * b + d - 2.0 / 3.0 * d)
t = ta - ma * x - ra / 2.0 * x ** 2 +
w1 / 24.0 * (x - a) ** 3 + w2 / 6.0 * (x - a - d) ** 3 + w3 / 24.0 * (x - a - 2.0 * b + d) ** 3 +
p1 * d ** 2 / 12.0 + p2 * (2.0 * b - d) ** 2 / 6.0 + p3 * d ** 2 / 12.0 +
p1 * d / 3.0 * (x - a - d) + p2 * (2.0 * b - d) / 2.0 * (x - a - 2.0 * b) + p3 * d / 3.0 * (x - a - 2.0 * b) +
p1 / 2.0 * (x - a - d) ** 2 + p2 / 2.0 * (x - a - 2.0 * b) ** 2 + p3 / 2.0 * (x - a - 2.0 * b) ** 2
v = va + ta * x - ma / 2.0 * x ** 2 - ra / 6.0 * x ** 3 +
w1 / 120.0 * (x - a) ** 4 + w2 / 24.0 * (x - a - d) ** 4 + w3 / 120.0 * (x - a - 2.0 * b + d) ** 4 +
p1 * d ** 2 / 12.0 * (x - a - d * 4.0 / 5.0) +
p2 * (2.0 * b - d) ** 2 / 6.0 * (x - a - d - (2.0 * b - d) * 3.0 / 4.0) +
p3 * d ** 2 / 12.0 * (x - a - 2.0 * b + d - d * 4.0 / 5.0) +
p1 * d / 3.0 * (x - a - d) ** 2 / 2.0 +
p2 * (2.0 * b - d) / 2.0 * (x - a - 2.0 * b) ** 2 / 2.0 +
p3 * d / 3.0 * (x - a - 2.0 * b) ** 2 / 2.0 +
p1 / 6.0 * (x - a - d) ** 3 + p2 / 6.0 * (x - a - 2.0 * b) ** 3 + p3 / 6.0 * (x - a - 2.0 * b) ** 3 + __g__ * (m - ma) * l ** 2
return q, m, t / ei, v / ei
end
def pt_mq5(tp, l, a, c, d, wf, t, h, e, g, κ, m1=0, m2=0, ka=0, kb=0, k1=0, k2=0, dn=10)
=begin
tp 支持条件
0 単純支持
1 Ａ端自由・Ｂ端固定
2 Ａ端固定・Ｂ端自由
3 両端固定
4 Ａ端ピン・Ｂ端固定
5 Ａ端固定・Ｂ端ピン
6 両端回転・鉛直バネ支持
l スパン m
wf 荷重 kN/m2
a 荷重の開始位置 m
c 荷重の終了位置 m
d 荷重の奥行 m
e ヤング率 kN/mm2
g 縦弾性係数 kN/mm2
κ せん断変形の形状係数
m1 Ａ端の回転外力 kNm
m2 Ｂ端の回転外力 kNm
ka Ａ端支持の回転バネ kNm/rad
kb Ｂ端支持の回転バネ kNm/rad
k1 Ａ端支持の鉛直バネ kN/mm
k2 Ｂ端支持の鉛直バネ kN/mm
dn 分割数
|d|
A _____wf B ___
_____／_____＼_____ _|_d
^ ^
|-a-|---2*b---|-c-|
0 L ---> x
|--------L--------|
q せん断力 kN
m 曲げモーメント kNm
t たわみ角 rad
v たわみ m
=end
ei = e * t * h ** 3 / 12.0 * 10 ** 6 # kN*m2
__g__ = (e / g) * (κ / 12) * (h / l) ** 2
ca, cb, qa, qb = cq5(wf, l, a, c, d, __g__)
ra, rb, ma, mb, ta, tb, va, vb =
boundary_condition(tp, l, ca, cb, qa, qb, ei, __g__, m1, m2, ka, kb, k1 * 1000, k2 * 1000)
pt = []
dx = l / dn
for i in 0..dn
x = dx * i
if x - dx < a && a < x && a != 0
q, m, t, v = mq5(a, l, a, c, d, wf, ra, ma, ta, va, ei, __g__, ka, kb, k1, k2)
pt << [a, q, m, t * 1000, v * 1000]
elsif x - dx < l - c && l - c < x && c != 0
q, m, t, v = mq5(l - c, l, a, c, d, wf, ra, ma, ta, va, ei, __g__, ka, kb, k1, k2)
pt << [l - c, q, m, t * 1000, v * 1000]
end
q, m, t, v = mq5(x, l, a, c, d, wf, ra, ma, ta, va, ei, __g__, ka, kb, k1, k2)
pt << [x, q, m, t * 1000, v * 1000]
end
return pt
end
tp, l, a, c, d, w, t, h, e, g, κ = ARGV.map { |x| eval(x).to_f }
pt = pt_mq5(tp, l, a, c, d, w, t, h, e, g, κ)
pt.each { |x| p ["x Q M T V", x[0], x[1].round(3), x[2].round(3), x[3].round(3), x[4].round(3)] }
__END__
画面に
と表示されます。
参考までに、ruby/tk なら
のように表示できます。
本稿のはりの計算式は小さなたわみに限られる実用式です。従来の計算式は曲げ変形にしか対応できませんでしたが、せん断変形の影響も考慮できます。書籍や論文で公開したものではありません。このブログではじめて紹介した計算式です。
よくよく考えると、当方は「はりの曲げ」を「フックの法則」に結び付けることを中学生のころから課題としてやっています。当時、平野くんが校庭で錆びた小さな帯鋸のかけらを拾い上げて、これに重りをつけて実験するから手伝ってくれと言われたことから始まっていて、おぼろげに学習発表会で平野くんが説明する横に立っていただけの記憶も残っています。
曲げ・せん断・剛域を考慮した「たわみ角法の基本式」は耐震設計ではよく知られてもいます。「たわみ角法の基本式」は「両端曲げのフックの法則」でもあり、ちからと変位が比例関係であれば、曲げであろうとせん断であろうとそれを加えたものであろうとたわみはたわみなので、オイラー梁とティモシェンコ梁が同じ計算式で解けても不思議はありません。なかったことが不思議なくらいです。
次回は 集中曲げを取り上げます。
