【ゲーム物理】離散要素法DFMで簡単物理エンジン！①　基礎編 #JavaScript

この記事は、基礎編 です。
この記事を読むことで、フルスクラッチでスイカゲームが作れるようになります。

高校物理を未履修でも実装はできる、という記述を心掛けます。
しかしベクトルの知識は求められます。（ゲームプログラマー向けベクトル学習記事はいずれ書きます）

はじめに

物理シミュレーション、自前の実装は難しいですよね。
高校物理の知識だけでは、運動方程式で例えば落下の軌道を再現できても（運動学）、例えば、スイカゲームのように力のつりあいを保ってボールが積み上げることは再現するのは難しいです。
その理由は、運動の自由度（用語）が高く、方程式の数をそこまで増やすのが難しいため、運動方程式を立式して連立方程式ソルバーへ入れ込んでも各接触力が一意に定まらないからです。（あと単純に立式をプログラミングするのがメンドウ）

その一方で、 記述が簡単 で、物理シミュレーション学域でも使用されているぐらい それっぽい 手法があります。
それは離散要素法(Discrete Element Method; DEM) です。

離散要素法とは？

学術的な詳細は、こちらを参照：https://www.jstage.jst.go.jp/article/jjspe/84/7/84_615/_pdf

超簡単に

分野の人に怒られるぐらい 超簡単 に書くと、

物体どうしの接触力の公式

があり、

それを運動方程式に代入する

だけです。
これだけで 意識しなくとも力のつり合いが保たれるシミュレーションができます。

簡単でしょ？

ちょっとだけ解像度を上げる

離散要素法では、 一瞬だけ物体同士がめりこむ ことを許します。

そしてめりこんだ瞬間に、 反発力（接触力）を作用させて、めりこみ状態から離れるように作用させます

イメージしやすい映像（接触力をゆるゆるにしています）

こんなゆるゆるの物理シミュレーションで大丈夫かよと心配される方のために、ちゃんとパラメーターを設定した映像を。

接触力の公式

事前に設定するべきパラメーター

ばね定数 $k$
めりこんだときに、どれくらいの強度で反発するか を調整するパラメーターです。
ずいぶん高めに設定されます。
反発係数（跳ね返り係数） $e$
衝突した直後に、直前と比べてどれくらいの速さで離れるか（正確には相対速度） を表します。
通常0~1の数値で、1だと、直前同じ速さで離れ、0.5だと半分、0だと止まります（ひっつきます）
質量 $m$
正確に言うと「力を受けたときの加速度の変化の比例係数」ですが、不正確に言うと「重さ」ですね。
運動量保存の法則により、 重い物体は衝突を受けても動きづらく、軽い物体はぶっとびやすい です。直観通りですね。これを念頭に設定してください。

計算するべき公式

接触力

F_{法線方向} = -k\delta -\eta v_{r, 法線方向}

です。

$F_{法線方向}$は、法線方向の接触力。
$k$ は ばね定数
先述の通り
$\delta$ は オーバーラップしている（めりこんでいる）度合い
僕は単純にめりこんでいる部分の長さ（下図の赤線の長さ）を代入しています。ちょっと複雑な図形だと体積の方がいいんだろうけど、求積の計算量が辛いかな。
$\eta$は 粘性係数
反発係数$e$から計算されます。　計算式は後述。
流体力学を心得ていないのであまりよく知らないです。
$v_{r, 法線方向}$は相対速度
$v_{r, 法線方向}= v_{相手, 法線方向} -v_{自分, 法線方向}$
あるいは
$v_{r, 法線方向} = (v_{相手}-v_{自分})の法線方向成分$
で求めます。
速度の法線方向成分を求めるには、正射影ベクトルを関数化しておくと便利です。

粘性係数

上記の粘性係数を求める公式は簡単。

\eta = -2\log{e}\sqrt{\frac{km}{(\log{e})^2 + \pi^2}}

で求められます。

$\eta$は粘性係数。
先述の通り。
$e$は反発係数。
先述の通り。ネイピア数とは違います！！
$\log{e}$は$e$の自然対数。
$e$はネイピア数ではないので$1$にはなりません。（限りません）
$k$はばね定数
先述の通り。
$m$は力を適用される物体の質量

接触力を算出したらどうするの？

運動方程式にそのまま入れたらOKです。

ここでは運動方程式をご存じでない方向けの説明です。

古典力学・ニュートン力学では、世の中の物体はこの方程式に従うとされています。

F=ma

$F$は力（のベクトル）
$m$は質量
$a$は加速度（のベクトル）

加速度$a$は速度の速度です
もし、速度が$3km/h$だと、「1時間後に3km移動している」という意味ですよね。
加速度が$3km/h^2$だと、「1時間後に速度が3km/hだけ速くなっている」という意味です。
簡単ですね。

気を付けるべきなのは$F$の計算。
受けている力ベクトル全てを合算しないといけません
たとえば物体1, 2 のみから力を受けていたら
$F = F_1 + F_2$
になります。

実装

接触力の算出

頭空っぽで、このような関数を書けばいいです。
JavaScriptの例です。

// 接触力を計算する関数
function CalculateContactForce(spring, restitution, overlap, speed_relative, mass){
    // 粘性係数を計算する関数
    function CalculateViscosity(spring, restitution, mass) {
        let temp = Math.sqrt(Math.log(restitution)**2 / (Math.log(restitution)**2 + Math.PI**2));
        //もし反発係数が（ほぼ）0だったら...
        if (Math.abs(restitution) < 1e-2) {
            //... INF/INF の不定形になるが、この場合は発散速度が同じなので1になる
            temp = 1
        }
        return 2 * temp * Math.sqrt(mass * spring);
    }

    const force_spring = spring * overlap;
    const vigosity = CalculateViscosity(spring, restitution, mass);
    const force_elastic = vigosity * speed_relative;

    return -force_elastic - force_spring;
}

運動方程式への適応

接触力は計算するだけでは意味がありません。運動方程式を通じて、実際に物体の運動に反映させる必要があります。

ここでは、動かす物体にこのようなクラスを付与していると仮定します。

class SimulatedObject{
    constructor(適切なもの) {
        this.position = 位置ベクトル
        this.velocity = 方向ベクトル
        this.mass = 質量//（数値）
        this.force_inframe = 方向ベクトル
    }
}

力を累積させる

運動方程式の説明で述べた通り、$F=ma$の$F$は合算された力ベクトルであるべきです。
そのため、フレームごとに受けた力全てを足し算しておき、最後に運動方程式を適用します。

すなわち、力を受けるたびに

//受けた力ベクトル
let force_got = /*ベクトルの大きさ*/CalculateContactForce(適切なパラメーター) * /*方向単位ベクトル*/normal_unit

simulatedObject.force_inframe += AddVec(simulatedObject.force_inframe, force_got)

という風に。AddVec()はベクトルの足し算の関数です。作ってください。

運動方程式の適用・移動

このように累積させた力を、運動方程式を通じて適用にします。

//運動方程式に基づいて加速度を計算
let acceleration = MultiplyVec(1/simulatedObject.mass, simulatedObject.force_inframe)

//force_inframeは用済みなので、次のフレームのためにリセット
simulatedObject.force_inframe = [0,0]

MultiplyVec()はスカラー値とベクトルの掛け算をする関数です。作ってください。

加速度を求めたら、このように物体を運動させます。

let delta_time = 前フレームと現フレームの時間

//加速度から速度を計算
simulatedObject.velocity = AddVec(simulatedObject.velocity, MultiplyVec(delta_time, acceleration))

//速度から座標を計算、移動
simulatedObject.position = AddVec(simulatedObject.position, MultiplyVec(delta_time, simulatedObject.velocity))

delta_timeの掛け算に注意してください。
加速度、速度は「1秒（単位時間）当たりの変化量」を表しますが、何も処理せずに足し算すると、暗黙に1秒経過したときの変化量を足すことになります。もしFPS（フレームレート）が60なら、現実時間1秒に、60秒分シミュレーションしてしまうことになります。

delta_timeの計算は、ストップウォッチを実装してフレームごとの経過時間を算出するのもよし、事前にFPSを設定している場合は逆数1/FPSのいわばSPFを代入してもよいです。

実装例

離散要素法を用いて、ニュートンのゆりかごを実装した例です。
レポジトリ：https://github.com/konbraphat51/NewtonsCradle
動作ページ: https://konbraphat51.github.io/NewtonsCradle/

ベクトルの計算がめんどくさければ、ぜひブラウザで動く弊ライブラリHotSoupScriptをご利用ください。

終わりに

できるようになったこと

ボール同士の力学シミュレーションを実装できるようになりました。

次回予告

四角形などに関しては、回転もあるため、今度はモーメントや角速度の計算が必要になります。

お願い

いいね頂けると泣きながら喜びます><

【ゲーム物理】離散要素法DFMで簡単物理エンジン！① 基礎編