Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

CA Tech LoungeAdvent Calendar 2023

Day 6
@konbraphat51inGoogle Developers Students Club: the University of Tokyo

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

Last updated at Posted at 2023-11-19

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

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

はじめに

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

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

離散要素法とは？

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

超簡単に

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

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

があり、

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

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

簡単でしょ？

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

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

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

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

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

接触力の公式

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

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

計算するべき公式

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

です。

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

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

\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をご利用ください。

終わりに

できるようになったこと

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

次回予告

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

お願い

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

3
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?