外部変形は データのやり取りをテキストファイルで行うので プログラム言語は 自由に選ぶことができます。図形は機能的かつシンプルなため、数多くのユーザーに受け入れられています。
maximaではりを解いています。
(はりを解く(両端回転ばね＋鉛直ばね：集中曲げ))
KB : 回転ばね Rotation Spring
KV : 鉛直ばね Virtical Spring
L : スパン長 Span Length
a : 集中曲げの位置 Position of load
x : はりの位置 Position
PM : 集中曲げ 左回りが正 turn-left is positive-direction
Q : せん断力 Shear ↓-↑ is positive-direction
M : 曲げモーメント 左回りが正 Bending Moment
T : たわみ角 Deflection Angle
V : たわみ Deflection
EI : 曲げ剛性 一定 Bending Stiffness (constant)
RA, MA, TA, VA : A 点の反力、曲げモーメント、たわみ角、たわみ
RB, MB, TB, VB : B 点の反力、曲げモーメント、たわみ角、たわみ
[ 実用弾性曲線式 ]
EI * V(4)(x) = -p(x)
p(x) = -w(x)
Q(x) = C1 + ∫ p(x) dx
M(x) = C2 + ∫ Q(x) dx
r(x) = -M(x) / EI
T(x) = C3 + ∫ r(x) dx
V(x) = C4 + ∫ T(x) dx
ただし
w(x) : 分布荷重 上方向が正
p(x) : 分布荷重 kN/m 上方向が正
Q(x) : せん断力 kN ↓-↑が正
M(x) : 曲げモーメント kNm 左回りが正
r(x) : 曲率 1/m 上に凸が正
T(x) : たわみ角 rad 左回りが正
V(x) : たわみ m 上方向が正
EI : 曲げ剛性(一定) kN*m*m
C1,C2,C3,C4 : 積分定数
スクリプト例
:はりを解く(両端回転ばね＋鉛直ばね：集中曲げ)
@echo off
path C:\maxima-5.49.0\bin;%path%
more +4 %0 | maxima --very-quiet & pause & goto:eof
/*
: PM
: _________/________
◎ ◎ KB
> > KV
: 0 a L ---> x
: |--------L--------|
y
|
|_____x
／
z
PM : 集中曲げ 左回りが正 kNm
a : 集中曲げの位置 m
KB : 回転ばね kNm/rad
KV : 鉛直ばね kN/m
Q : せん断力 ↓･↑ が正 kN
M : 曲げモーメント 左回りが正 kNm
T : たわみ角 rad
V : たわみ mm
L : スパン長 m
EI : 曲げ剛性(一定)
エンコード UTF-8(BOM無し) で保存して実行してください。
*/
assume(L>0,a>=0)$
P: 0$
define(Q1(x), C1)$
define(M1(x), C2 + integrate(Q1(x), x))$
define(T1(x), C3 - integrate(M1(x) / EI, x))$
define(V1(x), C4 + integrate(T1(x), x))$
define(Q2(x), C5)$
define(M2(x), C6 + integrate(Q2(x), x))$
define(T2(x), C7 - integrate(M2(x) / EI, x))$
define(V2(x), C8 + integrate(T2(x), x))$
co: [Q1(a) - Q2(a) = P, /* P は 上向きが正 */
M1(a) - M2(a) = -PM, /* PM は 左回りが正 */
T1(a) = T2(a),
V1(a) = V2(a)]$
s: solve(append(
[M1(0) =-T1(0)*KB, Q1(0) = V1(0)*KV, M2(L) = T2(L)*KB, Q2(L) =-V2(L)*KV], co),
[C1, C2, C3, C4, C5, C6, C7, C8])$
[Q1, M1, T1, V1]: ev([Q1(x), M1(x), T1(x), V1(x)], s)$
[Q2, M2, T2, V2]: ev([Q2(x), M2(x), T2(x), V2(x)], s)$
[C1,C2,C3,C4,C5,C6,C7,C8]:ev([C1,C2,C3,C4,C5,C6,C7,C8],s)$
Q:if x <= a then Q1 else Q2$
M:if x <= a then M1 else M2$
T:if x <= a then T1 else T2$
V:if x <= a then V1 else V2$
xMmax:a$
/*
RA: 6*KV*PM*(KB*(a-L)*a-EI*L)/((KB*L+6*EI)*KV*L^2+24*EI*KB)$
MA:-KB*PM*((KB*L*(a-L)*(3*a-L)+2*EI*(3*a*(a-2*L)+2*L^2))*KV*L-24*EI*(KB*(a-L)-EI)
)/((KB*L+2*EI)*((KB*L+6*EI)*KV*L^2+24*EI*KB))$
RB:-RA$
MB: KB*PM*((KB*L*a*(3*a-2*L)+2*EI*(3*a^2-L^2))*KV*L+24*EI*(KB*a+EI)
)/((KB*L+2*EI)*((KB*L+6*EI)*KV*L^2+24*EI*KB))$
TA:-MA/KB$
VA: RA/KV$
TB: MB/KB$
VB: RB/KV$
Q1: RA$
Q2:-RB$
M1: MA+RA*x$
M2: MB-RB*(x-L)$
T1: TA-MA*x/EI-RA*x^2/(2*EI)$
T2: TB-MB*(x-L)/EI+RB*(x-L)^2/(2*EI)$
V1: VA+TA*x-MA*x^2/(2*EI)-RA*x^3/(6*EI)$
V2: VB+TB*(x-L)-MB*(x-L)^2/(2*EI)+RB*(x-L)^3/(6*EI)$
*/
[L,a,PM,B,D,E]: [5,2,-50,0.15,0.30,6.5]$
KB:1000*10$
KV:1000*10$
EI: E * B * D^3 / 12 * 10^6$ /* kN*m*m */
/* xVmax 100,000,000 divisions */
vvv:0$ xxx:0$
m:L*100$ dx: L/m$
for i:0 thru m do
if abs(v:float(ev(V,x=dx*i))*1000) > abs(vvv) then [vvv, xxx]: [v, dx*i]$
xx1:xxx$ mxxx:xxx/dx$ m:m*100$ dx: L/m$ dt:1$ xxx:0$
for i:mxxx*(100-dt) thru mxxx*(100+dt) do
if abs(v:float(ev(V,x=dx*i))*1000) > abs(vvv) then [vvv, xxx]: [v, dx*i]$
xx2:if xxx=0 then xx1 else xxx$ mxxx:xxx/dx$ m:m*100$ dx: L/m$ dt:0.01$ xxx:0$
for i:mxxx*(100-dt) thru mxxx*(100+dt) do
if abs(v:float(ev(V,x=dx*i))*1000) > abs(vvv) then [vvv, xxx]: [v, dx*i]$
xx3:if xxx=0 then xx2 else xxx$ mxxx:xxx/dx$ m:m*100$ dx: L/m$ dt:0.0001$ xxx:0$
for i:mxxx*(100-dt) thru mxxx*(100+dt) do
if abs(v:float(ev(V,x=dx*i))*1000) > abs(vvv) then [vvv, xxx]: [v, dx*i]$
xVmax:if xxx=0 then xx3 else xxx$
/* 総合図 */
pm:if a > L-a then [[xMmax, 0], [xMmax, ev(M1,x=ev(xMmax))]] else [[xMmax, 0], [xMmax, ev(M2,x=ev(xMmax))]]$
pv:[[xVmax, 0], [xVmax, ev(V,x=ev(xVmax))*1000]]$
plot2d([Q,M,T*1000,V*1000,[discrete,pm],[discrete,pv]], [x,0,L], grid2d, [title, "Euler-Bernoulli's beam"],
[legend,"Q kN","M kNm","T x1/1,000 rad","V mm","Mmax","Vmax"],
[style,lines,lines,lines,lines,linespoints,linespoints],
[color,red,blue,green,magenta,blue,magenta])$
/* 画面表示 */
print("Virtical Deflection")$
m:L*100$ dx: L/m$
for i:0 thru m do
if mod(i,m/10) = 0 then printf(true, "V(~06,3f)=~09,3f~%",i*dx,float(ev(V,x=dx*i))*1000)$
print("Vmax of calcurated value")$
printf(true, "V(~09,06f)=~015,12f~%",ev(xVmax),ev(V,x=ev(xVmax))*1000)$
print("Bending Moment")$
m:L*100$ dx: L/m$
for i:0 thru m do
if mod(i,m/10) = 0 then printf(true, "M(~06,3f)=~09,3f~%",i*dx,float(ev(M,x=dx*i)))$
print("Mmax of calcurated value")$
printf(true, "1(~09,6f)=~015,12f~%",ev(xMmax),ev(M1,x=ev(xMmax)))$
printf(true, "2(~09,6f)=~015,12f~%",ev(xMmax),ev(M2,x=ev(xMmax)))$
printf(true, "A(~09,6f)=~015,12f~%",0,ev(M,x=0))$
printf(true, "B(~09,6f)=~015,12f~%",L,ev(M,x=L))$
?sleep(10)$
quit()$
(総合図)
KB=1000*10 kNm/rad
KV=1000*10 kN/m
L=5.0m
a=2.0m
PM=-50 kNm
B x D = 150x300mm
E65-F225 構造用集成材 すぎ(JAS)
(計算結果)
