Blenderで花火大会

2022/08/18

夏なのに引きこもりなので花火を作ります．

花火の作成

なにはともあれ，まずは花火の表現方法を考えます．結果だけ知りたい方はジオメトリノードを構築まで飛ばしてください．

採用しなかった方法：物理演算とモーションブラー

花火の運動を再現したい場合，真っ先に思いつくのはBlenderのパーティクル鋼体シミュレーションです．ちょっと実験してみましょう．

まぁこんな結果になるんですが，写真としては違和感があります．というのも，「花火大会」で画像検索してみると……

写真の花火は，「点」ではなく「線」で写っています．

これは，撮影時に露光時間を長く取っているためです．

開花した花火が消えるまでシャッターを開いておくと、（中略）花火の軌跡をきれいに捉えることができます。シャッター速度が短すぎると（中略）花火が点で写ってしまいます。

出典: 花火大会の彩りを撮る, ソニー

というわけで，「Blender　シャッタースピード」とか「Blender　長時間露光」で検索してみると，「モーションブラーを使え」との記述が出てきます．実際に使ってみると，

……あれ？なんだかガビガビです．

私の設定が悪いのかもしれませんが，どう頑張っても綺麗な結果は得られませんでした．

また，この方法では不便な点が多いように感じます．レンダリングしないと軌跡の形がわからないし，マテリアルの調整も大変そうです．

他にも，パーティクルの数を増やしたり，子パーティクルを設定したりして擬似的に軌跡を表現する方法も試しましたが，どれもきれいな結果にはなりませんでした……．

採用した方法：運動方程式＋ジオメトリノード

一番良いのは，鋼体の軌跡をメッシュとして直接得られることです．

そもそも，花火の動きは単純なので，時刻からその星（空中で光る火薬1つ）の位置を代数的に決定できるはずです．つまり，ステップバイステップの数値解析は必要ありません．

したがって，運動方程式を解いて星の位置を直接計算し，ジオメトリノードによって花火のメッシュを生成します．

運動方程式を解く

変数や関数を以下のように定義します．

$t$ ：花火が開く瞬間を $0$ とする時刻， $s$
$m$ ：星の質量， $kg$
$k$ ：空気抵抗係数，スカラー
加速度
- $\vec{a}(t)$ ：時刻 $t$ における星の加速度ベクトル， $(m/s^2, m/s^2, m/s^2)$
  - $a_x(t)$ ：時刻 $t$ における星の加速度ベクトルの $x$ 成分， $m/s^2$
- $\vec{g}$ ：重力加速度， $(m/s^2, m/s^2, m/s^2)$
速度
- $\vec{v}(t)$ ：時刻 $t$ における星の速度ベクトル， $(m/s, m/s, m/s)$
  - $v_x(t)$ ：時刻 $t$ における星の速度ベクトルの $x$ 成分， $m/s$
- $\vec{v}_0$ ：時刻 $0$ における星の速度ベクトル， $(m/s, m/s, m/s)$
  - $v_{0x}$ ：時刻 $0$ における星の速度ベクトルの $x$ 成分， $m/s$
位置
- $\vec{p}(t)$ ：時刻 $t$ における星の位置ベクトル， $(m, m, m)$
  - $p_x(t)$ ：時刻 $t$ における星の位置ベクトルの $x$ 成分， $m$
- $\vec{p}_0$ ：時刻 $0$ における星の位置ベクトル， $(m, m, m)$
  - $p_{0x}$ ：時刻 $0$ における星の位置ベクトルの $x$ 成分， $m$

星の運動方程式は重力と空気抵抗¹を用いて以下のようになります．

m \vec{a}(t) = m\vec{g} - k\vec{v}(t)

両辺を $m$ で割れば，

\vec{a}(t) = \vec{g} - \frac{k}{m}\vec{v}(t)

加速度 $\vec{a}(t)$ が得られます．

ここからは，ひとまず $x$ 成分に着目して解いていくことにします．

a_x(t) = g_x - \frac{k}{m} v_x(t)

簡単のために， $r = \frac{k}{m}$ とおきます．

a_x(t) = g_x - r v_x(t)

これを $t$ で積分して速度 $v_x(t)$ を求めます．変数分離を使います．

\begin{align*} a_x(t) &= g_x - r v_x(t) \\ \frac{d v_x(t)}{dt} &= g_x - r v_x(t) \\ \frac{d v_x(t)}{dt} &= -r \{v_x(t) - \frac{g_x}{r}\} \\ \frac{1}{v_x(t) - \frac{g_x}{r}} \frac{d v_x(t)}{dt} &= -r \\ \frac{1}{v_x(t) - \frac{g_x}{r}} d v_x(t) &= -r dt \\ \int \frac{1}{v_x(t) - \frac{g_x}{r}} d v_x(t) &= \int -r dt \\ \log \{v_x(t) - \frac{g_x}{r}\} &= -rt + C_1 \\ v_x(t) - \frac{g_x}{r} &= e^{-rt + C_1} \\ v_x(t) - \frac{g_x}{r} &= e^{-rt} e^{C_1} \\ v_x(t) &= e^{-rt} e^{C_1} + \frac{g_x}{r} \qquad (1.1) \\ \end{align*}

$C_1$ は積分定数です．

ここで， $v_x(0) = v_{0x}$ より

\begin{align*} v_x(0) &= v_{0x} \\ e^{C_1} + \frac{g_x}{r} &= v_{0x} \\ e^{C_1} &= v_{0x} - \frac{g_x}{r} \qquad (1.2) \\ \end{align*}

となるため， $(1.2)$ で $(1.1)$ を書き換えると

v_x(t) = e^{-rt} (v_{0x} - \frac{g_x}{r}) + \frac{g_x}{r}

となり，速度 $v_x(t)$ が得られました．

これを更に $t$ で積分することで，位置 $p_x(t)$ を求めます．

\begin{align*} v_x(t) &= e^{-rt} (v_{0x} - \frac{g_x}{r}) + \frac{g_x}{r} \\ \frac{d p_x(t)}{dt} &= e^{-rt} (v_{0x} - \frac{g_x}{r}) + \frac{g_x}{r} \\ p_x(t) &= \int \{e^{-rt} (v_{0x} - \frac{g_x}{r}) + \frac{g_x}{r}\} dt \\ p_x(t) &= \int e^{-rt} (v_{0x} - \frac{g_x}{r}) dt + \int \frac{g_x}{r} dt \\ p_x(t) &= (v_{0x} - \frac{g_x}{r}) \int e^{-rt} dt + \int \frac{g_x}{r} dt \\ p_x(t) &= -\frac{1}{r} (v_{0x} - \frac{g_x}{r}) e^{-rt} + \frac{g_x}{r} t + C_2 \qquad (2.1) \\ \end{align*}

$C_2$ は積分定数です．

ここで， $p_x(0) = p_{0x}$ より

\begin{align*} p_x(0) &= p_{0x} \\ -\frac{1}{r} (v_{0x} - \frac{g_x}{r}) + C_2 &= p_{0x} \\ C_2 &= \frac{1}{r} (v_{0x} - \frac{g_x}{r}) + p_{0x} \qquad (2.2) \\ \end{align*}

となるため， $(2.2)$ で $(2.1)$ を書き換えると

\begin{align*} p_x(t) &= -\frac{1}{r} (v_{0x} - \frac{g_x}{r}) e^{-rt} + \frac{g_x}{r} t + \frac{1}{r} (v_{0x} - \frac{g_x}{r}) + p_{0x} \\ p_x(t) &= \frac{1}{r} (v_{0x} - \frac{g_x}{r}) (1 - e^{-rt}) + \frac{g_x}{r} t + p_{0x} \\ \end{align*}

となり，位置 $p_x(t)$ が得られました．

$y$ 成分と $z$ 成分についても同様なので，時刻 $t$ における位置 $\vec{p}(t)$ は次のようになります．

\vec{p}(t) = \frac{1}{r} (\vec{v_0} - \frac{\vec{g}}{r}) (1 - e^{-rt}) + \frac{\vec{g}}{r}t + \vec{p_0} \\

ただし

r = \frac{k}{m}

です．

ジオメトリノードを構築

以下が最終的なジオメトリノードです．（画像をクリックすると拡大します）

このノードは，面から法線方向に星を射出します．そのため，球以外にも様々なジオメトリを入力することで，ほとんどの種類の花火を表現できます．

使用例

例えば次のような設定の花火ジオメトリノードを球に適用します．横向きの重力加速度は横風を意図しています．

ジオメトリノードの結果，以下のようなそれっぽいメッシュが得られます．

それっぽいですね．

解説

各部分について細かく説明していきます．

上図では，軌跡の辺を生成しています．「面にポイント配置」によって，星を配置する場所をランダムに決めています．「メッシュライン」によって，星の軌跡を表現します．なお，ポイントとメッシュラインそれぞれでIDを格納しておきます．

上図では，先程求めた $\vec{p}(t)$ を使って軌跡の位置を計算しています． $t$ はメッシュラインが生成した各頂点のIDを使って算出されます． $v_0$ は面のノーマルとノーマル速度から与えられます．また， $v_0$ を任意の係数でランダマイズする機能もあり，軌跡が単調になるのを防ぎます．