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

はりの実用式は
  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))
でした。

数式処理ソフト 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, κ]: [5, -10, 0.15, 0.3, 6.5, 6.5/15, 1.2]$      /* (M9) */
    [EI, GA]: [e*B*D^3/12, g*B*D] * 10^6$      /* (M10) */
    plot2d([Q, M, T*1000, V*1000], [x, 0, L], [legend, "Q", "M", "T", "V"])$      /* (M11) */
となりました。

ruby はフリーのプログラム言語です。スクリプト言語と呼ばれコンパイルは不要です。体感ですが、処理速度はエクセル並みです。

ruby はダウンロードしたのちインストールすれば使えます。

ruby はエディタで編集と実行ができるので 1980年代の BASIC に似た使い方ができます。エディタは実行機能を備えた Mery が便利です。

Mery で以下のスクリプトを書いて ファイル名 "000.bat" で保存して実行すれば図を描けます。支持条件は o---o、----|、|----、|----|、o---|、|---o に対応できます。
: 等分布荷重
@echo off
path C:\maxima-5.46.0\bin;%path%
ruby -x %~f0 4 5 -10 0.15 0.3 6.5 6.5/15 1.2
pause
goto:eof

#!ruby
def pt_qm(tp, l, w, b, h, e, g, κ, f1 = 'MS Gothic,8', f2 = 'MS Gothic,10', f3 = 'MS Gothic,10')
  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( if #{tp} = 0 then /* 単純支持 */
      [ev(M, x = 0) = 0, ev(V, x = 0) = 0, ev(M, x = L) = 0, ev(V, x = L) = 0]
    elseif #{tp} = 1 then /* A端自由・B端固定 */
      [ev(Q, x = 0) = 0, ev(M, x = 0) = 0, ev(T, x = L) = 0, ev(V, x = L) = 0]
    elseif #{tp} = 2 then /* A端固定・B端自由 */
      [ev(Q, x = L) = 0, ev(M, x = L) = 0, ev(T, x = 0) = 0, ev(V, x = 0) = 0]
    elseif #{tp} = 3 then /* 両端固定 */
      [ev(T, x = 0) = 0, ev(V, x = 0) = 0, ev(T, x = L) = 0, ev(V, x = L) = 0]
    elseif #{tp} = 4 then /* A端ピン・B端固定 */
      [ev(M, x = 0) = 0, ev(V, x = 0) = 0, ev(T, x = L) = 0, ev(V, x = L) = 0]
    elseif #{tp} = 5 then /* A端固定・B端ピン */
      [ev(T, 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) */
    set_plot_option(
    [title, "Timoshenko-Ehrenfest beam / Lite & Seen Lite"],
    [xlabel, "位置(x) m"], [ylabel, "たわみ(V) mm"], 
    [legend, "Q kN  ", "M kNm ", "T x1000 rad", "V mm  ", "p kN/m "], /* key 凡例 */
    [color, blue, red, green, magenta, black],
    [style, [lines,1], [lines,1], [lines,1], [lines,5], [lines,1]]
    )$ /* (M10-1) */
    txt: "
    set title font '#{f1}';
    set xlabel font '#{f2}';
    set ylabel font '#{f2}';
    set key font '#{f2}';
    set key top left;
    set label 'Rectangle' at 2.5,-5 font '#{f3}';
    set label sprintf('b  = %.3f m', #{b}) at 2.5,-7.5 font '#{f3}';
    set label sprintf('D  = %.3f m', #{h}) at 2.5,-10 font '#{f3}';
    set label sprintf('e  = %.3f N/mm2', #{e}) at 2.5,-12.5 font '#{f3}';
    set label sprintf('g  = %.3f N/mm2', #{g}) at 2.5,-15 font '#{f3}';"$ /* (M10-2) */
    plot2d([Q, M, T*1000, V*1000, p], [x, 0, L], [gnuplot_preamble, txt])$ /* (M11) */
    ?sleep(5)$ /* (M12) 5秒表示 */
    quit()$ /* (M13) 終了 */
  Maxima
  f.close
end
tp, l, w, b, h, e, g, κ = ARGV.map { |x| eval(x) } # m, kN/m, m, m, kN/mm2, kN/mm2
pt_qm(tp, l, w, b, h, e, g, κ)
__END__

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

計算結果を数値で確認できれば便利なのですが maxima でやるとかなり面倒な感じです。


等分布荷重で1端ピン・他端固定のたわみが最大値となる正確な位置(x)を求めようとすると案外大変です。
公式集にはせん断変形を考慮したたわみの計算式はないし、マトリクス変位法の応力解析ソフトで細かく自動分割するのもまぬけな話です。実際の設計ではラフな位置がわかればいいため、この問題にあまり意味はないのですが、当方はひっかかりました。血まみれになりましたが maxima はスパっと答えを出してくれました。

δ(x) = w*x*(x-L)*(2*(3*__g__+1)*x^2-(6*__g__+1)*L*x-(72*__g__^2+30*__g__+1)*L^2)/(48*EI*(3*__g__+1))

 

x = -(((27*__g__^3+27*__g__^2+9*__g__+1)^(1/3)*(1528823808*__g__^9+3630956544*__g__^8+3670769664*__g__^7+ 省略

 

たわみ(δ(x))の計算式があれば問題は解決しそうです。位置(x)は複素数の実部を取り出す

必要があるなど答えがでるだけのレベルでしたが数値を代入してみると答え自体は間違っていないようでした。


(M11)を以下のようにすれば DXF ファイルに変換できます。
    plot2d([Q, M, T*1000, V*1000], [x, 0, L], [gnuplot_preamble, txt], 
    [gnuplot_term, "dxf"], [gnuplot_out_file, "c:/jww/0000.dxf"])$ /* (M11) */
CAD に読み込んで計測すればなんとかなります。それでも、もっといい方法があるはずです。

 

次回は maxima で解いた結果を数値で確認するやり方を紹介します。