力積・反発・摩擦

目標

力積で衝突応答を計算し、反発係数と摩擦が何を決めるかを、数字と絵で確認します。このラボが終われば、物理エンジンの応答の層を、自分で作れるようになります。

なぜ重要なのか

硬い物体の接触時間は、1ミリ秒にもならないので、60分の1秒で回るシミュレーションでは、接触の過程を真似られません。そのため、物理エンジンは、過程を飛ばして、結果だけを計算します。衝突の前後の速度の関係を方程式に立てて、それを満たす速度の変化量を、一度に求めて適用します。

ここで、学んでおくべき慣習が2つあります。1つは、逆質量を保存することです。動かない物体を逆質量0にしておけば、特別な処理なしに、同じ式がそのまま回ります。もう1つは、離れつつあるなら、手を出さないことです。この検査を抜かすと、離れようとする物体をもう一度引き寄せて、表面に張り付きます。

ステップ

/root/impulseに、ツールボックスを置きます。
/root/impulse/impulse.pyに、力積の公式を作ります。
/root/impulse/out/03-conserve.txtに、運動量とエネルギーを書きます。
/root/impulse/out/bounce.csv、04-bounce.txt、bounce.pngに、跳ねる高さを書きます。
/root/impulse/out/05-friction.txtに、滑る箱の結果を書きます。
/root/impulse/out/06-tangent.txtに、接線方向の摩擦を書きます。
/root/impulse/out/balls.pngと07-balls.txtに、箱の中のボール5つの結果を書きます。

参考

絵は、nohup python3 -m http.server 8080 -d /root/impulse/out &で立てて、Webプレビューでhttp://localhost:8080/を開きます。
よくあるミスの1つ目は、相対速度がすでに正(離れつつある)なのに、力積を適用することです。物体が表面に張り付きます。
よくあるミスの2つ目は、摩擦で、速度を0で切らないことです。速度が負に下がって、箱が後ろへ進みます。
よくあるミスの3つ目は、接線方向を、法線の力積を適用したあとの速度で求めることです。値が変わります。

描画のツールボックスを置く

/root/impulse/gfxlib.pyを例のとおりに保存し、/root/impulse/check.pyでテストパターンを描いて、/root/impulse/out/00-check.pngを作ってください。パターンは、64x64の黒い背景に、(0,0)から(63,63)まで、白(255,255,255)の対角線を引き、そのあとに、(0,32)から(63,32)まで、赤(255,0,0)の横線を重ねて引いたものです。

このラボからは、PNGエンコーダーを作り直しません。最初のラボで手で作ったものと同じコードを、ツールとして配ります。ここで学ぶのは、ファイル形式ではないからです。

ラボのPodにはボリュームがないので、前のラボで作ったファイルが残っていません。そのため、ラボごとに、ツールボックスを置き直すことから始めます。

Canvas(w, h, bg)を作り、line(x0, y0, x1, y1, rgb)で2本の線を引いたあと、write_png(path)で保存してください。交差点(32,32)が赤になるには、横線をあとで引く必要があります。

このラボでは、跳ねる高さのグラフとボールの軌跡を描くのに使います。

力積の公式

/root/impulse/impulse.pyに、resolve(v1, inv_m1, v2, inv_m2, n, e)を作ってください。速度は2次元のタプル、nは1から2に向かう単位ベクトルです。相対速度の法線成分が、0以上なら何もせずそのまま返し、負なら、j = -(1+e)*v_rel/(inv_m1+inv_m2)を求めて、v1 -= j*inv_m1*n、v2 += j*inv_m2*nとした結果を、(v1, v2)として返してください。

v_rel = dot(v2 - v1, n)です。この値が正だということは、2つの物体がすでに離れつつあるという意味なので、手を出してはいけません。この検査を抜かすと、離れようとする物体をもう一度引き寄せて、表面に張り付かせます。

質量ではなく逆質量を受け取る理由は、動かない物体を、inv_m = 0として自然に表現できるからです。無限大を入れると、割り算のたびに特別な処理が必要です。

確認は、同じ質量の2つを、正面衝突させてみることです。e=1なら、速度が互いに入れ替わり、e=0なら、どちらも同じ速度になります。

何が保存され、何が失われるか

質量2の物体1(速度(3,0))と、質量1の物体2(速度(-1,0))が、法線(1,0)でぶつかる場合を、e=1、0.5、0について解いて、/root/impulse/out/03-conserve.txtに、p_before=、ke_before=、p_e1=、ke_e1=、p_e05=、ke_e05=、p_e0=、ke_e0=の8行を書いてください。運動量は、x成分の和、運動エネルギーは、0.5*m*v²の和です(小数点以下6桁)。

逆質量は、それぞれ0.5と1.0です。

運動量は、3つの場合とも同じでなければなりません。2つの物体が、大きさが同じで向きが反対の力積を受けるので、構造的にそうなります。値が異なって出たら、符号が間違っています。

運動エネルギーは、e=1のときだけ保存され、eが小さくなるほど減ります。失われたエネルギーは、現実では、音と熱と変形に行きます。

e=0なら、2つの物体の速度が同じになります。くっついて一緒に動くという意味です。

跳ねる高さはeの2乗で減る

ボールをy=5.0、v=0から落として、g=-9.8、dt=1/240で、2400ステップ回してください。1ステップは、v += g*dt、y += v*dtで、y < 0になったら、y = -y、v = -0.7*vとして跳ね返します。結果を、/root/impulse/out/bounce.csvに、step,y,vのヘッダーと2401行で書き、最初の3回の跳ね返りの最高の高さを、/root/impulse/out/04-bounce.txtに、peak1=、peak2=、peak3=として書いてください。そして、/root/impulse/out/bounce.png(256x256、黒い背景)に、px = int(step*255/2400)、py = 255 - int(y/5.0*255)で、色(120,200,255)の点を打ってください。

跳ね返りの最高の高さは、yが隣り合う2つの値より大きいか等しい地点(局所最大値)です。最初の3つを選べば済みます。

比率を計算してみてください。peak2/peak1とpeak3/peak2が、どちらもe² = 0.49の近くで出る必要があります。速度がe倍に減れば、上がる高さはeの2乗倍に減るからです。

グラフを見ると、山がだんだん低くなる、おなじみの形が出ます。この絵が、反発係数の意味を、1枚で見せてくれます。

y = -yで跳ね返すのは、めり込んだ分を反射させて返す方式です。0で切り詰めるよりも、エネルギーの損失が少ないです。

滑る箱はどこで止まるか

速度5.0で滑る箱に、クーロン摩擦(mu=0.3、g=9.8)をかけて、dt=1/240で回してください。1ステップは、v = max(0, v - mu*g*dt)、x += v*dtで、vが0になれば、止まったことになります。かかった時間と移動距離を、/root/impulse/out/05-friction.txtに、steps=、stop_time=、distance=の3行で書いてください。

摩擦は、速度に比例せず、一定の減速を与えます。そのため、速度が直線的に減って、時間と距離を、解析的に求められます。

t = v0 / (mu * g)        = 5 / 2.94 = 1.7007 s
d = v0² / (2 * mu * g)   = 25 / 5.88 = 4.2517 m

自分で回した値が、この式と1パーセント以内で合う必要があります。ちょうど同じにはならないのは、離散的なステップが、最後の断片を切り落とすからです。

max(0, ...)がないと、速度が負に下がって、箱が後ろへ進みます。摩擦が物体を押してやることになりますが、これは、実際のエンジンでもよくあるミスです。

摩擦も力積である

/root/impulse/tangent.pyに、bounce_with_friction(v, inv_m, n, e, mu)を作ってください。動かない表面(逆質量0)にぶつかる物体1つを扱います。法線の力積jn = -(1+e)*dot(v,n)/inv_mを求めて適用したあと、接線方向t = normalize(v - dot(v,n)*n)について、jt = -dot(v,t)/inv_mを求め、[-mu*jn, +mu*jn]で切り詰めた値を適用して、最終的な速度を返します。(v_out, jn, jt)の3つの値を返し、v=(3,-3)、inv_m=1、n=(0,1)、e=0.5、mu=0.4の場合を、/root/impulse/out/06-tangent.txtに、vout=vx,vy、jn=、jt=として書いてください。

接線方向は、速度から法線成分を引いた残りです。その残りが0なら(正面からぶつかった場合)、接線が定義されないので、摩擦を飛ばす必要があります。

jtを切り詰めることが、クーロン摩擦のすべてです。接線速度を完全になくすのに必要な値が、mu*jnを超えなければ、物体がその場で止まり(静止摩擦)、超えれば、切り詰めた値だけが適用されて滑ります(動摩擦)。条件文1つで、2つの摩擦がすべて出ます。

与えられた場合では、jn=4.5、切り詰める前のjt=-3、上限が0.4*4.5=1.8なので、jt=-1.8になります。最終的な速度は、(1.2, 1.5)です。

法線成分と接線成分は、jnを適用する前の速度で計算します。順序を混ぜると、値が変わります。

箱の中のボール5つ

半径0.5のボール5つを、1辺が10の箱の中で、g=-9.8、dt=1/240、壁の反発係数0.9で、1800ステップ回してください。初期状態は、順番に、位置が(2,8)・(5,9)・(8,7)・(3,5)・(7,4)、速度が(3,0)・(-2,1)・(0,0)・(4,2)・(-3,-1)です。10ステップごとに、各ボールの位置を、/root/impulse/out/balls.png(256x256、黒い背景、sx=int(x*25.6)、sy=int(255-y*25.6))に、色(255,100,100)・(100,255,100)・(100,150,255)・(255,220,80)・(220,120,255)で打ち、/root/impulse/out/07-balls.txtに、b1=x,y,vx,vyからb5=までと、energy_start=、energy_end=を書いてください。エネルギーは、sum(0.5*(vx²+vy²) + 9.8*y)です。

壁の処理は、軸ごとに別々に行います。x - r < 0なら、x = r、vx = -e*vxで、反対側も同じです。yも同じです。

1ステップの順序は、重力の適用、位置の更新、壁の処理です。順序を変えると、ボールが壁を突き抜けます。

エネルギーは、必ず減る必要があります。反発係数が1より小さいので、壁にぶつかるたびに、少しずつ失います。増えたなら、壁の処理で位置を戻す部分が抜けているか、符号が間違っています。

軌跡が絵に残ります。放物線の断片がつながって跳ねる様子が、目に見えます。