[力学]

ウェブページ用の HTML言語の制約によって、このページの数式は、 のぞましくない形で表現されているところがある。

• 分数を斜線でしめしている (横線のほうがよい)。
• 平方根を sqrt(...) であらわす (√ をつかうべき。)

斜面をころがる円板の運動の計算

2021-12-23 に [剛体の力学 (2)]のページの中にふくめた 「小課題」の部分を独立ページにし、2022-01-06 に 補足する。

例、斜面をころがる円板 (この部分「剛体の力学 (2)」と重複)

質量 M、半径 a の円板が、 斜面 (かたむき θ) の最大傾斜方向に、すべらずにころがり落ちるとする。

• [(*) 2-18図 斜面を転落する円板]

円板にはたらく力は、重心に(地球の)重力 M g、 斜面との接点で 摩擦力 F と 垂直抗力 R。

x軸を斜面にそってくだる方向にとる。(水平方向でないことに注意。)

重心の運動の式 (2.53) は つぎのようになる。

M d²X/dt² = M g sinθ - F ... (a)

M d²Y/dt² = R - M g cosθ ... (b)

回転運動の式 (2.54) は つぎのようになる。

I_G d²φ/dt² = - F a ... (c)

F と R は「束縛力」であり、ここでは、 円板が斜面にめりこむこともすべることもなくころがるという条件 (Y = 0, dX/dt = - a dφ/dt)からきまる。

Y = 0 なので、(b) から、R = M g cosθ。

「慣性モーメントの計算」の「例2」から、I_G = (1/2) M a²。 これを (c) に代入すると、F = - (1/2) M a d²φ/dt²。

ところが、dX/dt = - a dφ/dt なので、a d²φ/dt² = - d²X/dt²。

したがって、F = (1/2) M d²X/dt²。

これを (a) に代入すると

M d²X/dt² = M g sinθ - (1/2) M d²X/dt²

したがって

d²X/dt² = (2/3) g sinθ

dX/dt = (2/3) g t sinθ

X = (1/3) g t² sinθ

cf. 摩擦がなく、ころがらずにすべるばあいは、d²X/dt² = g sinθ

X = X₁ まで行ったとする (教科書の l (エル小文字) を X₁ と書いておく。)

円板の高さは、初期の位置よりも X₁ sin θ だけ低くなっている。

円板の位置エネルギーは、初期よりも M g X₁ sinθ だけへっている。

時刻は t = sqrt((3 X₁/(g sinθ)) である。

重心の運動 (剛体の並進運動) の運動エネルギーは

(1/2) M (dX/dt)² = (2/3) M g X₁ sinθ ... (d)

重心のまわりの回転の運動エネルギーは (dX/dt = -a dφ/dt なので)

(1/2) I_G (dφ/dt)² = (1/2) (M a²/2) (dφ/dt)² = (1/4) M (dX/dt)²

とくに X = X₁のとき

(1/2) I_G (dφ/dt)² = (1/3) M g X₁ sinθ ... (e)

このばあい、位置エネルギー、並進運動のエネルギー、回転運動のエネルギーの合計が保存している。

一様な円板や円筒がすべらずにころがるときは、 並進運動と回転運動の運動エネルギーの比は、2:1 である。 この比は、剛体の形によってちがう。

小課題

斜面をころがる円板の例題で、M, a の値を適当にあたえて、 時間 t を等間隔にサンプリングして、 各時刻の、X, φ、水平・鉛直の座標軸からみた重心の位置 (X cosθ, -X sinθ)、 円周上につけた印の位置 (X cosθ + a cosφ, -X sinθ + a sinφ) を表にしてみよう。 できれば、水平座標を横、鉛直座標を縦にとった平面上で、重心と印の軌跡を図にしてみよう (Excelなどの表計算ソフトウェアのグラフのうち「散布図」の機能でできるはず。)。 g は便宜上 10 m/s² としてもよい。

注意

運動方程式の数値計算では、微分を有限時間間隔の差分におきかえて近似解 (いわゆる「数値解」)をもとめることがよくある。 今回の課題はそうではなく、微分方程式のまま積分によってもとめた解 (いわゆる「解析解」) の数値を見る計算である。 プログラムあるいはワークシートでは、 解析解の式に、有限の (無限小ではない) 時間間隔で数値を代入した結果を、 有限のけた数で近似して表示している。 解析解に数値を代入する計算式は、同じ時間ステップの他の数量の値を参照することがあるが、 ひとつ前の時間ステップの数量の値を参照することはない。 (数値解をもとめる式は、前のステップの値を参照するのがふつうである。)

表計算でやる例

ワークシート1 (「小課題」の範囲の作業)

Excel の「セルに名まえをつける」機能をつかったほうが計算式が見やすくなるはずだが、 まだためしていない。

ワークシートの左側、定数を置く

ワークシートの右側、時間変化する量を計算する

「↓」のところは、上のセルから、計算式をコピーする (数値をコピーするのではない)。 コピーすると、計算式のなかの、絶対参照(「$」つき)のセル番号はかわらないが、 相対参照(「$」なし)のセル番号は ずれていく。

時間ステップ数を、たとえば 50 までにするならば、「↓」のところは、ワークシートの55行めまでコピーをくりかえす。

各列ごとに、「セルの書式設定」の「表示形式」を「数値」にして「小数点以下の桁数」を 設定 (たとえば「3」に) したほうがよい。

円周上の印の軌跡のグラフを書いてみる

(操作例は Windows上の Excel 2013 のばあい)

ワークシート1のうち、M～N列、5～55行めを選択する。

「挿入」から「グラフ」、そのうちから「散布図」、 そのうちから「散布図 (直線とマーカー)」を選択する。

ワークシート2: 速度、角速度、エネルギーを追加

各時刻の 重心の速度 dX/dt、角速度 dφ/dt、 並進の運動エネルギー、回転の運動エネルギー、位置エネルギーを計算する (12月23日に「発展」として書いておいた項目)

ワークシート1をコピーしてワークシート2とする。

ワークシート1で計算した列のうちには、ここで不要なものもあるが、 便宜上、ワークシート1のA～N列はそのままにして、その右に列を追加する形にした。 必要のない列の「列幅」をせまく設定すれば、見たい複数の列をいっしょに見ることはできる。

列幅のつごうでつぎのような省略をしている。 「重心の速度」は重心の速度の斜面下向き成分 dX/dt である。 「E」は「エネルギー」である。

時間ステップ数を、たとえば 50 までにするならば、「↓」のところは、ワークシートの55行めまでコピーをくりかえす。

各列ごとに、「セルの書式設定」の「表示形式」を「数値」にして「小数点以下の桁数」を 設定 (たとえば「3」に) したほうがよい。

プログラム言語 Python でやる例

12月23日の小課題の範囲の作業 (「ワークシート1」と同様な計算)

時間変化する量の列の数値 (各時刻ごとに1行) を書きだす。

(現在のバージョンのプログラムでは、 各列の数量が何であるか、その単位は何かを、文字で書きだすことはまだできていません。)

• korogaru.py (文字コード UTF-8)

操作のごく概略

Windows に Anaconda を入れたばあい、Anaconda Power Shell のプロンプトから (「cd」で korogaru.py を置いたディレクトリを作業ディレクトリとして指定したあと)

python korogaru.py > korogaru.txt

を実行すれば、結果が korogaru.txt にはいる。 Linux でも同様だが、python3 というコマンドがはいっているばあいは、 python でなく python3 をつかうほうがよい。

korogaru.py のソースプログラム

(どこで何をやっているかを講義で説明する予定)

# korogaru.py
# 2021-12-23 増田 耕一
# ころがる円柱の運動
# 重心の位置と、円周上の印の位置
# 暫定版プログラム (各列の量と単位の説明を書きだしていない)

import math

# 定数を設定する
grav     = 10 # 重力加速度 [m/s2]
thetadeg = 30 # 斜面の傾き [度]
theta    = thetadeg * math.pi / 180  # 斜面の傾き [rad]
costheta = math.cos(theta)
sintheta = math.sin(theta)
mass     = 10 # 円板の質量 [kg]
radius   =  1 # 円板の半径 [m]
deltat   =  0.05 # 出力の時間間隔 [s]
nt       = 100   # 出力ステップ数

# 時間ステップ計算
for kt in range(nt+1):
    tim        = kt*deltat                       # 時刻 [s]
    xdownslope = grav * (tim**2) * sintheta / 3  # 重心の位置 (斜面下向き) [m]
    phi        = -xdownslope / radius            # 回転角 [rad]
    direction  = phi*180/math.pi % 360           # 円周上の印の方向 [度]
    hpos_center =  xdownslope * costheta         # 重心の水平位置 [m]
    vpos_center = -xdownslope * sintheta         # 重心の鉛直位置 [m]
    hpos_mark   = hpos_center + radius * math.cos(phi)  # 印の水平位置 [m]
    vpos_mark   = vpos_center + radius * math.sin(phi)  # 印の鉛直位置 [m]
    print(f'{kt:5}'+f'{tim:10.3f}'+f'{phi:10.3f}'+f'{direction:10.3f}'+\
          f'{hpos_center:10.3f}'+f'{vpos_center:10.3f}'+\
          f'{hpos_mark:10.3f}'+f'{vpos_mark:10.3f}')

円周上の点の軌跡の作図

• korogaru_plot1.py (文字コード UTF-8)

日本語文字のあつかいは、コンピュータのOSなどに依存するところがあるので、 ひとまず、表示させる文字列は英語にした (日本語ローマ字でもよいのだが)。

ソースプログラムの korogaru.py からの変更点

import math
import matplotlib.pyplot as plt
.....
# 定数を設定する
.....
# 配列変数を準備する
hpos_mark_list = []
vpos_mark_list = []

# 時間ステップ計算
for kt in range(nt+1):
.....
    hpos_mark   = hpos_center + radius * math.cos(phi)  # 印の水平位置 [m]
    hpos_mark_list.append(hpos_mark)
    vpos_mark   = vpos_center + radius * math.sin(phi)  # 印の鉛直位置 [m]
    vpos_mark_list.append(vpos_mark)
.....
    (printの実行文 (継続行をふくむ) は頭に「#」を打ってコメントにしておく)

# 作図
plt.plot(hpos_mark_list, vpos_mark_list)
plt.title('Trajectory of a point on the perimeter of a rolling disk')
# (ころがる円板の円周上の点の軌跡)
plt.xlabel('Horizontal position [m]')
# (水平位置)
plt.ylabel('Vertical position [m]')
# (鉛直位置)
# 画面に表示 (つぎの plt.show を生かす)
## plt.show()
# ファイルに保存 (つぎの4行を生かす)
imagefilename = "korogaru_plot1.png"
plt.savefig(imagefilename)
plt.clf()
plt.close()

結果の図

速度、角速度、エネルギーを追加 (「ワークシート2」と同様な計算)

各時刻の 重心の速度 dX/dt、角速度 dφ/dt、 並進の運動エネルギー、回転の運動エネルギー、位置エネルギーを計算する (12月23日に「発展」として書いておいた項目)

• korogaru2.py (文字コード UTF-8)

ソースプログラムの korogaru.py からの変更点

# 定数を設定する
.....
m_inert  = mass * (radius**2) / 2 # 慣性モーメント [kg m^2]

# 時間ステップ計算
for kt in range(nt+1):
.....
    vdownslope  = grav * tim * sintheta * (2/3)   # 重心の速度 (斜面下向き) [m/s]
    ang_speed   = -vdownslope / radius            # 角速度 [rad/s]
    e_transl    = mass * (vdownslope**2) / 2      # 並進の運動エネルギー [J]
    e_rotat     = m_inert * (ang_speed**2) / 2    # 回転の運動エネルギー [J]
    e_pot       = mass * grav * vpos_center       # 位置エネルギー [J]
    e_total     = e_transl + e_rotat + e_pot      # 力学的エネルギー [J]
    print(f'{kt:5}'+f'{tim:10.3f}'+f'{phi:10.3f}'+f'{direction:10.3f}'+\
          f'{hpos_center:10.3f}'+f'{vpos_center:10.3f}'+\
          f'{hpos_mark:10.3f}'+f'{vpos_mark:10.3f}'+\
          f'{vdownslope:10.3f}'+f'{ang_speed:10.3f}'+\
          f'{e_transl:10.3f}'+f'{e_rotat:10.3f}'+\
          f'{e_pot:10.3f}'+f'{e_total:10.3f}')

発展

(12月23日に「発展」としたうち、まだプログラム例をしめしてないことがら)

回転せずにすべり落ちるばあい (摩擦はないとする) と比較してみよう。