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

:::note
この記事は、**基礎編** です。
この記事を読むことで、フルスクラッチで[スイカゲーム](https://www.youtube.com/watch?v=zBjjg44Shmo&t=708s)が作れるようになります。

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

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

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

## 離散要素法とは？
:::note
学術的な詳細は、こちらを参照：https://www.jstage.jst.go.jp/article/jjspe/84/7/84_615/_pdf
:::

### 超簡単に
分野の人に怒られるぐらい **超簡単** に書くと、
* **物体どうしの接触力の公式**

があり、
* **それを運動方程式に代入する**

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

簡単でしょ？

### ちょっとだけ解像度を上げる
離散要素法では、 **一瞬だけ物体同士がめりこむ** ことを許します。

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

イメージしやすい映像（接触力をゆるゆるにしています）
![tama_20231119_202247.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602537/cb132b68-413d-0ea8-5046-9f9e41eb79c3.gif)

こんなゆるゆるの物理シミュレーションで大丈夫かよと心配される方のために、ちゃんとパラメーターを設定した映像を。
![tama_20231119_202708.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602537/67410c6c-a424-2553-2689-4f97c131aea4.gif)

### 接触力の公式
#### 事前に設定するべきパラメーター
* [ばね定数](https://wakariyasui.sakura.ne.jp/p/mech/dannsei/hukku.html) $k$
    めりこんだときに、**どれくらいの強度で反発するか** を調整するパラメーターです。
    ずいぶん高めに設定されます。
* [反発係数（跳ね返り係数）](https://wakariyasui.sakura.ne.jp/p/mech/hannpatu/keisuu.html) $e$
  衝突した直後に、直前と比べて**どれくらいの速さで離れるか（正確には相対速度）** を表します。
  通常0~1の数値で、1だと、直前同じ速さで離れ、0.5だと半分、0だと止まります（ひっつきます）
* [質量](https://www2.nhk.or.jp/kokokoza/watch/?das_id=D0022150204_00000) $m$
    正確に言うと「力を受けたときの加速度の変化の比例係数」ですが、不正確に言うと「重さ」ですね。
    [運動量保存の法則](https://www.try-it.jp/chapters-8001/sections-8222/lessons-8223/)により、 **重い物体は衝突を受けても動きづらく、軽い物体はぶっとびやすい** です。直観通りですね。これを念頭に設定してください。
    
#### 計算するべき公式
##### 接触力
```math
F_{法線方向} = -k\delta -\eta v_{r, 法線方向}
```
です。
* $F_{法線方向}$は、**法線方向の接触力**。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602537/ec1b44f1-4090-a9a4-1eef-acf37a9496cb.png)
* $k$ は **ばね定数** 
先述の通り
* $\delta$ は **オーバーラップしている（めりこんでいる）度合い**
僕は単純にめりこんでいる部分の長さ（下図の赤線の長さ）を代入しています。ちょっと複雑な図形だと体積の方がいいんだろうけど、求積の計算量が辛いかな。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602537/8515b7db-ea1e-083b-0c72-91ae6a2e96be.png)
* $\eta$は **粘性係数** 
    **反発係数$e$から計算されます。**　計算式は後述。
    流体力学を心得ていないのであまりよく知らないです。
* $v_{r, 法線方向}$は[相対速度](https://www.try-it.jp/chapters-8001/sections-8066/lessons-8073/)
    $v_{r, 法線方向}= v_{相手, 法線方向} -v_{自分, 法線方向}$ 
    あるいは
    $v_{r, 法線方向} = (v_{相手}-v_{自分})の法線方向成分$ 
    で求めます。
    速度の法線方向成分を求めるには、[正射影ベクトル](https://manabitimes.jp/math/933)を関数化しておくと便利です。
    ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2602537/d2218455-b131-2738-44ed-9def5975ca0a.png)

##### 粘性係数
上記の**粘性係数**を求める公式は簡単。
```math
\eta = -2\log{e}\sqrt{\frac{km}{(\log{e})^2 + \pi^2}}
```
で求められます。
* $\eta$は粘性係数。
    先述の通り。
* $e$は反発係数。
    先述の通り。**ネイピア数とは違います！！**
* $\log{e}$は$e$の自然対数。
    $e$はネイピア数ではないので$1$にはなりません。（限りません）
* $k$はばね定数
    先述の通り。
* $m$は力を適用される物体の**質量**

### 接触力を算出したらどうするの？
運動方程式にそのまま入れたらOKです。

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

古典力学・ニュートン力学では、世の中の物体はこの方程式に従うとされています。
```math
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`の例です。
```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;
}
```

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

ここでは、動かす物体にこのようなクラスを付与していると仮定します。
```javascript
class SimulatedObject{
    constructor(適切なもの) {
        this.position = 位置ベクトル
        this.velocity = 方向ベクトル
        this.mass = 質量//（数値）
        this.force_inframe = 方向ベクトル
    }
}
```

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

すなわち、**力を受けるたびに**
```javascript
//受けた力ベクトル
let force_got = /*ベクトルの大きさ*/CalculateContactForce(適切なパラメーター) * /*方向単位ベクトル*/normal_unit

simulatedObject.force_inframe += AddVec(simulatedObject.force_inframe, force_got)
```
という風に。`AddVec()`はベクトルの足し算の関数です。作ってください。

#### 運動方程式の適用・移動
このように累積させた力を、運動方程式を通じて適用にします。
```javascript
//運動方程式に基づいて加速度を計算
let acceleration = MultiplyVec(1/simulatedObject.mass, simulatedObject.force_inframe)

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

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

加速度を求めたら、このように物体を運動させます。
```javascript
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](https://github.com/konbraphat51/HotSoupScript)をご利用ください。

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

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

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