― 62 ―煙のビジュアルシミュレーション法に関する検討
籾山静香＊・竹下大樹

Study on Visual Simulation of Smoke
Shizuka MOMIYAMA ＊ and Daiki TAKESHITA（平成23年11月25日受理）
In this paper, we made and ran the two-dimensional fluid simulator for visual simulation ofsmoke. It is necessary to solve the simultaneous equations in a part of the processing. The iterativemethod is used as a method of numerical calculation. Jacobi method, Gauss-Seidel method, and SOR（Successive Over-Relaxation）method was adopted for solving simultaneous equations. We triedto specify the method which is the most eﬃcient to simulate by comparing the number of repetition
to convergence. As a result, SOR method is the best one because it takes the least calculation cost.

1.

緒言

自然界に存在する水や煙，雲など様々な流体の現象の CG による表現は流体のシミュレーションを基本として行うことができる。しかし流体のシミュレーションは計算量が多いのが難点である。
近年，インタラクティブなコンテンツの生成手法に関する注目が高まっており，流体シミュレーションをコンテンツに利用するにあたり，処理の効率化はますます重要となっている。
本研究では文献［1］に示される手法を元に 2 次元の流体シミュレーションプログラムを作成し，シミュレーションを行った。このシミュレーションの処理の一部では連立方程式を解く必要があり，数値計算の手法として，反復法を用いている。本研究では反復法からヤコビ法，ガウス・ザイデル法，SOR法を採用してシミュレーションを行い，収束までの反復回数を比較することで，流体のシミュレーションを効率的に行うことが可能な手法の特定を試みた。
2.

煙のビジュアルシミュレーション法

本手法では 2 次元空間において流体の運動シミュレーションを行い，煙密度と温度を移流することで煙の表現を試みる。シミュレーションでは流体を非＊秋田高専専攻科学生圧縮性流体として扱う。2 次元空間を格子で区切り，格子の中心に速度，圧力，温度，煙密度を定義し，各格子点における時間毎の値を更新する。速度と圧力の更新には，次の式で表される非圧縮性流体の運動方程式であるナビエ・ストークス方程式を用いる。
∂u＝－（u・∇）u－∇p＋v∇2 u＋F∂t（1）∇・u ＝0（2）ここで，v は動粘性率，F は外力を表す。また，式（1）の右辺は第一項から順に移流項，圧力項，粘性項，外力項をそれぞれ表している。実際のプログラムでは外力項，移流項，粘性項，圧力項の順に計算を行い，速度場と圧力場を更新する。外力には渦度閉じ込め，重力，浮力を考慮した。その後得られた速度場を用いて，温度の移流，煙密度の移流を行う。温度の移流と煙密度の移流の方程式は以下を用いた。
∂T＝－（u・∇）T
∂t（3）

∂ρ
＝－（u・∇）ρ∂t（4）ここで，ρは煙密度，T は温度を表す。
2.1 外力項
本研究では，外力として，渦度閉じ込め，重力，
平成24年２月

― 63 ―
煙のビジュアルシミュレーション法に関する検討浮力を考慮した。渦度閉じ込めとは流体に渦の消失を抑えるような力を加えて渦の維持を助けるもの
で，式（5）で表すことができる。
f ＝εh（N×W）

粘性項の計算は安定性を重視し，陰解法を用いる。
この方程式を陰解法によって，t から Δt で差分化す
ると次式のようになる。

（5）

ここでεは渦度閉じ込めの係数，h は格子幅であり，N，W は以降の式であらわされる。
ηN＝｜η｜（6）η ＝∇｜W｜（7）

W＝∇×u（8）図1

これをオイラー法によって時間刻みごとに速度を求める。
また，重力と浮力は式（9）によって表される。
fbuay＝－αρz＋β（T－Tamb）z（9）ここで，α および β は正の定数，z は z ＝（0，0，1）
で鉛直上向きのベクトルをそれぞれ現している。この式の第一項は重力によって下向きに働く力を示し，第二項は温度 T と環境温度 Tamb の差によって生じる浮力であり，上向きに働く力を示している。
2.2 移流項
移流は流体の持つ速度によって，着目する流体塊における量が輸送される現象である。この量とは速度，温度や煙密度を対象とすることができる。本研
究では 2 次元空間を格子で区切っているため，移流を求めるには，各格子点の量を更新することになる。
移流項の計算ではラグランジアン法のアプローチを用いて移流の計算を行う。これは流体を流体塊として扱い，格子点xの速度u（x, t）を用いて，流れを遡った位置 x' ＝ x－u（x, t）Δt での量を格子点 x での新しい値としてコピーする方法である。量 q を更新するには次の式を使う。
q（x, t＋Δt）＝ q（x－u（x, t）Δt, t）（10）

位置 x' は格子点上に一致するとは限らないため，位置 x' に最も近い 4 つの格子点の量を用いて線形補間し，結果を対象とする格子点 の値として更新する。
2.3 粘性項
流体は流れに対して速度の拡散が生じる。これを粘性といい，粘性項は∂u/∂t＝v∇2 u で表される。
秋田高専研究紀要第47号

流体の移流計算u（x, t＋Δt）－ u（x, t）＝ v∇2 u（x, t＋Δt）
Δt（11）

この式を u（x, t＋Δt）で整理すると，（I－vΔt∇2 ）u（x, t＋Δt）＝ u（x, t）（12）となる。ここで I は単位行列を表している。この式を中心差分で差分化する。ここで，u（x, t）＝bu（x, t＋Δt）＝u（13）（14）とし，式（12）に空間差分を用いると以下の式となる。
b＝u－vΔt  u i＋1, j －2u i, j＋u i－1, j u i, j＋1－2u i, j＋u i, j－1＋（Δx）2（Δy）2（15）ここでΔx ＝Δy として式を整理すると，式（16）が導き出される。

u i, j ＝（Δx）2u i＋1, j＋u i－1, j ＋ u i, j＋1＋u i, j－1＋ vΔt b（Δx）2vΔt ＋4（16）式（16）をヤコビ法によって，十分な回数計算することで，粘性項の計算を行うことができる。ヤコビ法の反復計算の回数を k として，式（16）において b＝u i, j とすると，式（17）のようになる。

― 64 ―
籾山静香・竹下大樹

u（k＋1）＝i, j（Δx）2（k）
（k）
（k）
u（k）
u i, ji＋1, j ＋ u i－1, j ＋ u i, j＋1＋ u i, j－1＋ vΔt（Δx）2vΔt ＋4p（k＋1）＝i, j（k）
（k）
（k）
p（k）
i＋1, j ＋ p i－1, j ＋ p i, j＋1＋ p i, j－1＋4－Δx2∇・w
Δt（17）（22）となる。この式（17）に 3 つの反復法を適用しシミュレーションを行う。となる。この式（22）に 3 つの反復法を適用しシミュレーションを行う。

2.4 圧力項
本研究では非圧縮性流体を仮定した。その場合，流体は非発散の条件∇・u ＝0 を満たす必要がある。
外力，移流，粘性の計算によって得られた速度場をw とすると，この速度場は非発散の条件を満たしていない状態となっているので，速度場を調整する必要がある。文献［2］によると，この速度場の調整は以下のように行われる。

3.

u＝w－Δt∇p（18）

反復法

3.1 ヤコビ法
ヤコビ法について，圧力項を例に説明する。ヤコビ法のプログラムでは図 2 のように，計算に利用する p（k）の値を格納する配列（領域 1）と，計算結果の p（k＋1） の値を格納する配列（領域 2）を 2 つ用意して計算を行うことになる。p（k＋1）
を求める場合にi, j（k）
は領域 1 の斜線部の p の値を式（22）に代入し計算を行い，結果を領域 2 の対応した領域に格納する。

速度場 u は圧力項を計算した後に得られる速度場であり，非発散の条件を満たす。速度場 w は外力，移流，粘性によって得られた速度場である。この式によって，速度場を更新するためには圧力を求める必要がある。圧力を求める式は以下となる。
1∇2p＝ ∇・wΔt（19）となる。これは圧力のポワソン方程式である。式（19）を解いて，圧力を導出し，式（18）によって，速度場を更新する。
この式の左辺を中心差分によって差分化すると，
pi＋1, j －2pi, j ＋ pi－1, j pi, j＋1－2pi, j ＋ pi, j －1 1＋＝ ∇・w
2Δt（Δx）2（Δy）（20）となる。ここでΔx ＝Δy として整理すると式（21）が求められる。

図2

ヤコビ法

3.2 ガウス ･ ザイデル法
ガウス・ザイデル法の場合には既に求められている値はただちに使用するという方法であるので，値を格納する領域は 1 つで済む。図 3 のように領域の左上からプログラムが走査していくものとすると，p（k＋1）の値が既に計算されている場合はその値を使用する。図 3 の太枠内の p（k＋1）
を求める場合を例にi, j（Δx）2∇・w
pi＋1, j ＋ pi－1, j ＋ pi, j＋1＋ pi, j －1－Δt
pi, j ＝
4（21）
式（21）をヤコビ法によって，十分な回数計算することで，圧力を導出することができる。ヤコビ法の反復計算の回数を k とすると，式（21）は

図3

ガウス・ザイデル法

平成24年２月

― 65 ―
煙のビジュアルシミュレーション法に関する検討

とると，4 近傍のうち，2 つの値は既に更新されている p（k＋1）の値を使用することになる。斜線部の値を式（22）に代入し計算を行い，結果を同じ領域に格納する。
3.3 SOR 法
SOR 法はガウス ･ ザイデル法を元にした手法である。手順を図 4 に沿って説明する。

この手順を領域 1 の全ての p（k）に対して繰り返す。
3.3.1 加速係数
SOR 法で使用するωとは，加速係数のことを指す。この加速係数とは一度の反復計算における計算処理量を加速する係数のことである。その式に適した加速係数はそれぞれ異なり，大きければ良いわけでも，小さければ良いわけでもない。適した加速係数を設定することにより計算手順を減少させることが出来る。
本研究では，1.0から2.0まで0.1刻みで加速係数を変化させ，収束までの反復回数を比較し，流体のプログラムに最も適した加速係数の選定を試みた。結果は表 1 に示すとおりである。
表1ω

加速係数選定結果
圧力

粘性合計

平均

1.016375.457

1.116235.410

1.215485.160

1.31597

5.323

1.416135.377

1.51734

5.780

1.61906

6.353

1.72255

7.517

1.8311710.3091.9617520.583

合計

平均3001.00

表 1 からわかるように，ω＝1.2のときに最も反復回数が少なくなっている。よって本研究ではω＝1.2と定めた。

図4SOR 法

①まず，領域 1 から領域内の全ての p（k）を領域 2 にコピーする。
②領域 1 の求めるべき格子点の圧力を p（k）
i, j として，ガウス・ザイデル法で計算する。このとき，4 近傍の値に既に SOR で計算された p（k＋1） の値があれば，使用する。
③②で求めたガウス・ザイデル法を用いた p（k＋1）のi, j値にωを掛けたものと，①でコピーした領域 2 のp（k）
i, j に（1－ω）を掛けたものを足し合わせた値が SOR 法を用いて計算した p（k＋1）となる。
i, j
秋田高専研究紀要第47号

3.4 収束判定条件
これらの反復法では収束判定条件を定める必要がある。収束判定条件とは，x（k）の値が真の解に十分に近づいたかを判定する条件のことである。これが満たされたとき反復計算を途中で打ち切る。収束判
定条件にはいくつかの方法があるが，本研究では以下のものを採用した。
ωをある小さな正数であるとする。そして x（k）に関する漸化式を k＝0, 1, 2…の順に計算していき，NN｜x（k＋1）－x（k）｜x（k＋1）｜iii ｜<εΣΣ
i＝1i＝1（23）となったときに計算を打ち切り，x（k＋1）を答えとする。実際の計算においては，x は圧力 p や速度 u が対応する。

― 66 ―
籾山静香・竹下大樹

4.

描画方法

今回流体のシミュレーションで用いた描画方法について以下の図を用いて説明する。
（1） を煙密度の値を持つ格子点，太線を境界とする。

以上の手順を踏むことで得られた輝度値を用い，OpenGL を用いて，煙の描画を行う。図 5 の左図のような煙密度を設定した場合，図 5 の右図のような描写となる。

（2）周囲 4 点の格子点 の値の平均を求め，その値を用いて輝度値を点 に設定する。

図55.

（3）周囲 4 点のする。の輝度値の平均を求め， に設定

結果

流体のシミュレーションのプログラムにおいて，粘性項と圧力項に反復法の 3 つの手法を適用した
際の反復回数の結果を表 2 に示す。また100，150，200ステップ目のシミュレーション結果を図 5 に示す。
表2
手法をつなぎ合わせ，それぞれ三角形反復回数の平均

ヤコビ法

ガウス・ザイデル法

SOR 法

圧力6.573

5.457

5.160

粘性1.0

1.01.0

計算対象（4）点 と点を生成する。

描写の例

表 2 を見ると，粘性項の計算では空気の粘性係数が非常に小さく反復計算を必要としないため，手法による差は現れなかったが，圧力項の計算ではSOR 法が最も反復回数が少なく収束しているのがわかる。また，図 6 のシミュレーション結果を見ると，煙の広がり具合と煙密度に差が出ているのが分かる。以下でそれらについての検討を述べる。
平成24年２月

― 67 ―
煙のビジュアルシミュレーション法に関する検討

表3
手法
計算対象

反復回数の平均ヤコビ法

ガウス・ザイデル法

SOR 法

10e－2

6.573

5.457

5.160

10e－67197.78

4623.49

1867.18は非常に多くの反復計算が行われる。これではインタラクティブなコンテンツを作成するツールの開発を目的とした場合，計算コストが高い。よって，本研究では反復回数を減らすことを優先的に考え，ε＝10－2 としてシミュレーションを行った。
次に煙密度に差が出ていることについては，移流の計算手法が影響しているものと考えられる。図 8は150ステップ目のシミュレーション結果と圧力分布図を並べたものである。

図66.

シミュレーション結果（150ステップ目）

考察

まず，煙の広がり具合に差が出ていることについては，圧力項の計算結果の違いが流体の速度に影響しているものと考えられる。今回シミュレーションに用いた収束判定条件の式（23）ではε＝10－2 を採用しているため，完全に計算が収束する前に計算を打ち切っている。実際ε＝10－6 を採用し計算を行った場合，結果は図 7 のようになり，手法による差はほとんど見られない。

図7

図8

シミュレーション結果と圧力分布図

圧力分布図の中央にある黒い領域に注目すると，ヤコビ法では広い領域で緩やかに黒色に変化しているのに対して，ガウス・ザイデル法や SOR 法では，狭い領域で急激に黒色へ変化している。これは領域が狭く色が濃いほど，強い圧力勾配が発生していることを表している。強い圧力勾配が発生していると，黒色領域の外側から内側へ大きな速度が生じる。
2.2. で説明したように，今回のシミュレーションでは移流の計算手法により，流れを遡った位置の値を用いて計算を行う。よって煙のない場所の値を用いて煙のある位置の計算を行ったために，煙密度が薄くなってしまったものと考えられる。

収束判定条件ε＝10－6（150ステップ）7.
しかし，ε＝10－6 とした場合，行われる反復計算は表 3 のようになる。
表を見るとわかるように，ε＝10－6 とした場合に秋田高専研究紀要第47号

結論

本研究では煙のビジュアルシミュレーション法の効率化を図るために，2 次元の流体シミュレーショ

― 68 ―
籾山静香・竹下大樹

ンプログラムを作成し，いくつかの数値計算の手法を適用，比較することで，効率的にシミュレーションを行うことのできる手法の検討を行った。その結果以下のことがわかった。

今後の課題としては，CG 法（共役勾配法）やICCG 法（不完全コレスキー分解付共役勾配法）などを適用してみることで，より最適な手法を選定することや，計算手法を変えることによって得られる画像の違いについて検討を行うことなどが挙げられる。

（1）本研究の 2 次元の流体シミュレーションプロ
グラムにおいて，SOR 法を用いて圧力を求めるには，ω＝1.2が最も適した加速係数である。 8. 参考文献
（2）ヤコビ法，ガウス・ザイデル法，SOR 法で比較すると，SOR 法が最も反復回数が少なく収 ［1］J. Stam. Stable Fluid. In Proceedings of ACM
SIGGRAPH 99, Computer Graphics Proceedings,
束する。
Annual Conference Series, 121-128, 1999.
以上のことから流体のシミュレーションに用いる ［2］R. Fedkiw, J. Stam and H. Jensen. Visual Simula数値計算の手法は，3 手法の中で SOR 法が計算コス
tion of Smoke. In Proceedings of ACM SIGGRAPH 2001, Computer Graphics Proceedings,
トの点で最も優れているといえる。
収束判定条件における係数ε＝10－2 とε＝10－6 のAnnual Conference Series, 15-22, 2001.
［3］Hubert Nguyen :「GPU Gems3」，ボーンデジ場合のシミュレーション結果を比較すると，煙の形
タルに差が出ているが，本研究ではインタラクティブな［4］髙橋大輔：「数値計算」，岩波書店
コンテンツの作成を目的としているため，計算効率を重視し，現状では許容範囲であると考えている。
［5］越塚誠一：「数値流体力学」，培風館

平成24年２月