はりの計算式を解く(4)

はりの実用式は
　　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)
　　V' (x) = T(x) = C3 + ∫ r(x) dx　　　　　　　　　　　　 (5)
　　V (x) = V(x) = C4 + ∫ T(x) dx + κ * EI / GA * ∫ Q(x) dx　 (6)
　　ただし
　　κ　　　せん断変形の形状係数
　　EI　　　曲げ剛性(一定)
　　GA　　せん断剛性(一定)
　　x　　　位置
　　L　　　スパン長
　　p(x)　　分布荷重
　　w　　　等分布荷重
　　Q(x)　　せん断力
　　M(x)　　曲げモーメント
　　r(x)　　曲率(δ''=φ=r(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))
でした。

Mery で以下のスクリプトを書いてファイル名 "000-1.bat" で保存して実行すれば図が描けることも紹介しました。

: 等分布荷重単純支持
@echo off
path C:\maxima-5.46.0\bin;%path%
ruby -x %~f0 5 -10 0.15 0.3 6.5 6.5/15 1.2
pause
goto:eof

#!ruby
def pt_qm_tp_0(l, w, b, h, e, g, κ)
f = open "| maxima --very-quiet", "w"
f.print <<~Maxima
p: -w$ /* (M1) */
Q: C1 + integrate(p, x)$ /* (M2) */
M: C2 + integrate(Q, x)$ /* (M3) */
r: -M$ /* (M4) */
T: C3 + integrate(r, x)$ /* (M5) */
V: C4 + integrate(T, x) + κ * EI / GA * integrate(Q, x)$/*(M6)*/
s: solve([ev(M, x = 0) = 0, ev(V, x = 0) = 0, ev(M, x = L) = 0, ev(V, x = L) = 0]
, [C1, C2, C3, C4])$ /*(M7)*/
[Q, M, T, V]: subst(s, [Q, M, T/EI, V/EI])$ /* (M8) */
[L, w, B, D, e, g, κ]: [#{l}, #{w}, #{b}, #{h}, #{e}, #{g}, #{κ}]$ /* (M9) */
[EI, GA]: [e*B*D^3/12, g*B*D] * 10^6$ /* (M10) */
plot2d([Q, M, T*1000, V*1000, p], [x, 0, L])$ /* (M11) */
?sleep(5)$ /* (M12) ５秒表示 */
quit()$ /* (M13) 終了 */
Maxima
f.close
end
l, w, b, h, e, g, κ = ARGV.map { |x| eval(x) } # m, kN/m, m, m, kN/mm2, kN/mm2
pt_qm_tp_0(l, w, b, h, e, g, κ)
__END__

(M1)～(M13)が maxima です。

今回は、はりを１０等分した位置のたわみの数値を調べようと思います。

maxima で (12)の後に以下２行を追加します。
pv: makelist([x, ev(V*1000)], x, 0, L, L/10)$ /* (M12-1) */
write_data(float(pv), "C:/jww/ztemp_pv.txt")$ /* (M12-2) */
テキストファイル "C:/jww/ztemp_pv.txt" に位置とたわみのデータを保存しました。

ruby で
File.foreach("C:/jww/ztemp_pv.txt") do |line|
p ["x, V(x)", *line.split.map { |x| x.to_f }]
end
とするとテキストファイル "C:/jww/ztemp_pv.txt" をコマンドプロンプトに