プログラミングと数値計算, 物理学

質点と壁との反発を表す運動方程式

コメントする

質点が壁に衝突し、反発することを数式で表現します。

  1. 壁の定義
  2. 運動方程式の導出
  3. プログラム
  4. 実行結果
  5. 参考文献

壁の定義

壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出

まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式

ポテンシャル\(V(x)\)保存力下の質点の運動を考えます。
壁が
\(f(x,t)=0\)
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数\(c\)を用いて、

と書くことが出来ます。
ここで、時刻\(t＝t_i\)は、
\(
f(x(t),t)=0~~\to~~t=t_i
\)
を満たす時刻(壁との衝突時刻)です。\(f(x,t)\)の時間微分は、

と書くことが出来ます。
\(c\)を定めるために式(2)の両辺を、時刻\(t=t_j\)周りを微小時間\(\Delta\)で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻\(t\)で成立するとします。計算を進めると、

という結果が導かれます。ここで、\(\delta_{i,j}\)はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、\(c\)を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度\(v(t_i+0)\)が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度\(v(t_i-0), v(t_i+0)\)と、壁の衝突前後の速度\(v_\text{w}(t_i-0), v_\text{w}(t_i+0)\)の間には、反発係数\(e\)を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、\(e\)を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置\(x(t_i)\)における壁の速度\(v_\text{w}(t_i)\)を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置\(x\)と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻\(t=t_i\)の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式

壁

と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻\(t=t_j\)周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、\(e\)は反発係数、\(\mathbf{v}_\text{w}(t)\)は壁の速度ベクトル、\(\mathbf{n}\)は壁の法線ベクトルであり、\(\mathbf{n}\)は

と表されます。
実際に数値計算を行う際には、衝突時刻\(t_i\)と位置\(\mathbf{r}=\mathbf{r}_i\)を求めた後、

に従って速度ベクトルを変更すると良いと思います（根拠は特にありません）。

多次元の反発


衝突時刻\(t=t_i\), 点\(\mathbf{r}(t_i)\)で質点が\(f(\mathbf{r},t)=0\)で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ\(\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),\)と置きます。
位置\(\mathbf{r}_i\)における壁の単位法線ベクトルを\(\mathbf{n}, (\mathbf{n}^2=1)\)と書くと、

と書くことが出来ます。ここで、\(c\)は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、\(e\)は反発係数で\(0\leq e\leq 1\)である(\(e=0\):完全非弾性衝突、\(e=1\):完全弾性衝突。
また、\(v_\perp, v_{\text{w}\perp}\)は\(\mathbf{v}, \mathbf{v}_\text{w}\)の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の\(c\)を求めましょう。式(27)の両辺に\(\mathbf{n}^\mathsf{T}\)を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が\(f(\mathbf{r},t)=0\)を満たす線と表現されていれば、時刻\(t=t_i\)で\(\mathbf{n}\)は

と求める事が出来ます。
続いて壁面の速度\(\mathbf{v}_\text{w}(t)\)を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の\(\mathbf{r}(t)\)は質点の位置ベクトルです。すなわち、衝突が起こる時刻\(t=t_i\)の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム

2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
\(
f(x,y,t)=\pm g(x,y,t)
\)
は同じです。

プログラム

gnuplot用のスクリプト

実行結果

プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

\(
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
\)

重力下の運動、動かない壁


\(
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
\)

壁
\(
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
\)

重力下の運動、動く壁


\(
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
\)

参考文献

3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方

関連記事

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です