セルオートマトンを用いた砂時計の CG シミュレーションに関する研究
2004MT032 礒部 浩誠

2004MT063 水野 志保
指導教員

1 はじめに
現在, 進歩した CG の技術を利用し, 規模の大きい自然物や物理現象が再現されてきている. その中でも水や空気などの流体を用いたシミュレーションは多く研究されているが, コンピュータの発展により, 最近では流体を無数の粒子の集合として扱うことができるようになった.
流体に比べて粉体の研究は少ないが, 技術の発展によりコンピュータは大きな計算能力をもつようになった. これにより, 様々な装置内で運動する粒子挙動が見られるようになり, 粉体現象の研究や論文などに有効なツールとして粉体シミュレーションが注目されている. 砂時計
は砂粒の崩壊, 落下, 蓄積などの様子を観察できるのでシミュレーションの題材として非常に興味深い. しかし, 砂時計のシミュレーションを行っている研究はさほど多くないので, 我々はこれに取り組むことにした. 本研究では砂時計における砂粒の落下と減少, それに伴う崩壊や蓄積といった砂の挙動をシミュレーションすることを目的とし, セルオートマトンを用いて, 砂の蓄積, 崩壊を二次元で表現する数値モデルを作成した.

2004MT075 岡田 佳子

金 知俊
可能であるものとし, 対象とする粉体を多数の粒子を含む粉体要素の集合で表現する. 次に, それぞれの粉体要素の運動をラグランジュ的に追跡し,粉体挙動のシミュレーションを行なう方法.
4. セルオートマトン法
格子状のセルと単純な規則からなる離散的計算モ
デル. このモデルは非常に単純化されたモデルであるが複雑な自然現象を表現することができる.

n+1nj-1

jj+1

有限差分法

有限要素法

SPH法

セルオートマトン

2 粉体の計算手法
粉体のシミュレーションを行なう上で考えられる手法
は有限差分法, 有限要素法, SPH 法, そして, 本研究で用いたセルオートマトン法がある. 以下にこれらの手法について説明する (図 1).
2.1 計算手法の比較
1. 有限差分法 (Finite Diﬀerence Method, FDM)
微分を差分商で近似することにより微分方程式を解く方法. 流体の基礎方程式で記述されている一次導関数をテイラー展開から導き出された式によって表すことで離散化を行う手法.
2. 有限要素法 (Finite Element Method, FEM)
複雑な形状, 性質を持った連続体を要素で分割し計算する方法. 領域全体の支配方程式を小領域に分割し, 各小領域には比較的単純で共通な補間関数を用いながら全体の補間精度を上げることが可能である. 解析領域全体の弱形式は積分で表されるため, 各々の要素内の積分の総和としても表せる. 利点は, 複雑な形状であってもメッシュに分割できれば解析ができるという点である.
3. SPH 法 (Smoothed particle hydrodynamics)
流体を粒子要素で分割し, これらの要素と運動が連続体力学に基づく支配方程式を満たすようにして処理する. 流動する粉体は連続体的取り扱いが

図1

様々な計算手法

我々は他の手法に比べ単純な手順であり, 計算時間も少なくてすむことからセルオートマトンを選んだ.
2.2 セルオートマトン法
同じ大きさの正方形 (セル) で区切られた面があり, 各セルに白または黒の石をおくことを想定する. 隣のセルの状態によって, 各セルの状態が変わるとしてその変わり方の規則を定める. 各時間において石の並べ方を決め,その操作を繰り返すとある時間からこれ以上変わらないパターンとなる. この手順がセルオートマトン法の原型であり, その特徴は以下のようになる.

• 同じ大きさの均一なセルを想定する• 白と黒の 2 値だけではなく, 各セルは k 種類の状態 (色) をとることができる• 次の時間のセルの状態は現在の状態と隣り合うセルの状態だけの局所的な規則による• でき上がるパターンは初期のセル状態と適用する規則によって決まる
セルオートマトン法は非常に単純な手順と簡単な規則から成り立っているが, 複雑なパターンを生み出すことができる. セルを横一列に並べたものが一次元モデルで,一次元セルの並び方の初期状態と状態変化の規則を設定し, 各セルの時間ステップ後の変化の例を図 2 に示す. これにより得られた結果を縦に並べると模様が現れ, そのパターンは二次元の特徴を備えている.

• 隣り合うセルのみを比較• L が既定の高さ以上のとき崩壊• 崩壊する場合 L/2 が隣のセルに移動

ステップ013315421(a) 砂山オートマトン

2

図2L/2

一次元セルオートマトン法によるセル配列の時間変化例
Lまた, 一次元のセルオートマトンモデルを一般化すると, セルの位置座標を i, 時間を t とし, 各セルの状態は 0から k − 1 までの k 個のうちから一つをとる. 位置 i のセルの時間ステップ t での状態を ait で表す. セル i の次の時間ステップ t + 1 の状態 ait+1 は, 隣りのセルとして,ai−rから ai+rまでの 2r + 1 個を考慮し, 隣りのセルのtt
状態を変数とする関数 F によって与えられる.
i−r+1ait+1 = F (ai−r, · · · , ait , · · · , ai+r−1, ai+rt , attt ) (1)式 (1) セル i の時間ステップ t + 1 の状態は近傍のセルの時間ステップ t の状態から決定される [2].

3 砂時計のモデル化
本研究では, 砂時計の挙動を以下の 3 領域に分けてモデル化し, シミュレーションを行う.

• 砂粒の減少による流れ• 砂粒の落下• 砂粒の堆積による崩れ

iijj(b) 崩壊の様子
図 3 一次元砂山オートマトン

本研究では, 砂時計の上部, 下部それぞれを一次元砂山オートマトンで表現する. セルオートマトンを用いることによって, 連続体力学に基づいた粒子間の相互作用の力などは考慮されなくなるが, セルに含まれる粒子の動きは追うことができる. 砂時計の上部で考えてみると,1粒の砂が落ちると実際にはその空いたスペースに様々な場所から砂粒が入り込み複雑な力の働きや挙動をする.
それによって膨大な計算時間が必要になる. それに対し, セルオートマトンでは領域が格子で分かれているので, 砂粒が落下した場合, その列の領域の砂の量が少なくなるだけなので, 容易に砂の挙動を求めることができる.
従って, 他の粉体シミュレーションの手法に比べて計算時間が少なくてすむ (図 4).

3.1 砂粒の流れと崩壊
3.1.1 粒子の挙動
本研究では, 砂の挙動のシミュレーションによる砂時計の再現を目的としている. そのため, 砂粒の流れと崩壊の挙動を表すモデルとして, 一次元砂山オートマトンを用いて粉体の複雑な運動パターンを表現する. 一次元砂山オートマトンは単純な粉体オートマトンモデルの代表的なもので, 一次元のセルに砂柱の高さを状態量として格納し, 隣り合う砂柱の高さの差 L が限界値を越えると,安定を失って左右どちらかに崩れ, 隣のセルとの高さの差を小さくするように砂粒の再配分が行われるルールを設定したものである (図 3)[3]. セルオートマトンを用いれば, 単純なルールで複雑な粉体の現象を再現することができる. 本研究でも文献 [1] の手法を元に単純なルールで砂時計の砂を再現することにした.
本研究で砂時計を表すために用いている一次元砂山
オートマトンのルールは, 以下の 3 点である.

一般的な粒子の移動、落下

図4

セルオートマトンによる粒子の移動、落下オートマトンの利点

境界線
砂時計の容器の範囲内でのみ砂粒の移動を行うようにするため, 容器を表す境界線を設定し, 上部の砂柱の長さが境界線を上回っている場合と下部の砂柱の長さが境界線を下回っている場合にセル間での砂粒の移動を行う.
なお, 砂時計の容器の境界線を描く際にはシグモイド曲線を利用した (図 5).

3.1.2

y=a1 + be−kx(2)式 (2) を逆関数にすると以下の式が得られる.

y=1a−xlog+hKbx(3)よって式 (3) より, 境界線を設定した. 但し, a, b, k, hは任意の係数および定数である.
y=ay(1) = a * k / 4y(2) = 0

逆関数にして砂時計型に(x0, a/2)
b * e -k * x0 = 1

粒が上部から下部に向けて流れ落ちる速度は, 容器内に残っている砂粒の量に無関係に一定である [4]. これに基づき, 一定時間毎に一定量の砂粒を上部から落下させた.

4 シミュレーション結果
3 章のモデルをプログラムとして実装し, シミュレーションを行った. 画像表示には EGGX ライブラリを用いた [5]. シミュレーションの結果を図 7 に示す.

y=0

a = 1, b = 1, k = 2

図5

シグモイド関数

また, 境界線と砂柱が重なることにより, セルに格納されている値が整数ではなくなる. これより, セルに格納されている値は砂粒の個数ではなく, 砂柱の高さということになる.
3.2 砂粒の落下
砂時計上部からの砂粒の落下は自由落下運動を用いる. 砂粒は自由落下により下降していき, 下の砂山の頂上にぶつかった時点で砂山に融合する. この砂粒と砂山の衝突を以下のステップにより判定する (図 6).

• ステップ 1
砂粒が上部から落下したらそのセルから落下した分の高さを引く• ステップ 2
砂粒の位置 y を自由落下運動の式により求める• ステップ 3
その時点での下部の砂山の中心部の高さ d を確認• ステップ 4y ≥ d ならば衝突しない → ステップ 2• ステップ 5y < d ならば衝突する → 下部の砂山の中心部の高さに堆積した分を足す

落下した分減少(a) 0秒後

図7(c) 180秒後シミュレーション結果

計算条件
本研究では, 高さ 10cm, 横幅 5cm の砂時計を想定している. 砂柱のセル配分は 101 分割で, その中で中央の 3つのセルから砂が落下するようになっている. また, セルの高さの変化量は横幅の 1/3 ずつとし, L が崩壊する値が 2 であることから, 安息角 (砂山が崩壊しない最大傾斜角度) は 18.43◦ と求められる. 時間刻みは 0.05 秒とする. なお, 砂の総量とは, 各セルに格納されている高さの合計である. それを擬似的に砂の総量とみなす.
4.2 形状・傾斜・時間の変化による挙動の比較
本研究では, 容器の形状として (a) シグモイド関数を用いた砂時計, (b) 傾斜角 45.00◦ の直線を用いた砂時計,(c) 二次曲線を用いた砂時計の 3 種類の 3 分時計を用意した (図 8). この 3 種の砂時計の動きを比較してみると,ほぼ同じような動きをした. また, それぞれの砂の総量を比較した結果を表 1 に示す. これより, 砂の落ちきる時間は容器の形状に左右されず, 砂の総量で決定されることが確認できた.

4.1表1

砂粒の位置y(b) 90秒後

形状の違いによる砂粒の総量の比較

形状(a) シグモイド(b) 直線(c) 二次曲線

砂の総量

14094.83

14238.00

14077.60

砂山の位置d

図6

砂粒の落下の様子

なお, 上部から砂粒が落下する際, 砂時計の容器の側壁にかかる圧力は容器内の砂粒の量によらないので, 砂次に, 傾斜角による速度の違いを観察するために, 一番
傾斜角の違いが分かりやすい (b) で, 表 1 と砂の量を同じ値に設定し, 容器の傾斜角が 18.43◦ , 33.69◦ の場合を比較した. それぞれの砂が落下しきる時間を比較した結果を以下の表 2 に示す. これより, 容器の傾斜角がゆるやかになると砂粒が落下しきる時間が長くなることが分は実際の砂時計に近い挙動を再現できている. しかし, 2は, 実際の砂時計の挙動とは異なる動きである. 3.1.1 節
では, 図 4 のような砂粒の流出の仕方の違いによる計算時間の削減という利点を述べているが, その違いによって 2 のような実際の挙動と異なった結果が生まれてしまった.

5 おわりに(a)シグモイド関数(b)直線

図8(c)二次曲線

3 種類の形状による比較

かる. これは, 砂柱の崩壊の仕方による違いである. 3.1.1
節で述べたように, 実際の砂時計と本研究のシミュレー
ションには図 4 のように砂粒の流出の仕方に違いがある. これより, 実際の砂時計では砂粒が一定の間隔で流れ続けるのに対し, 本研究で行ったシミュレーションでは砂山の勾配が安息角に到達した部分で一時的に崩壊が停止してしまい, それによって落下部分への砂粒の供給が不足し落下が断続的になる場合があり, 流出速度が低下してしまったからと考えられる.
表2傾斜角の違いによる落下時間の比較

傾斜角18.43◦

33.69◦

45.00◦

落下時間276.70

200.85180.25

さらに, 時間の違いを観察するために, (a) の砂時計の1 分時計と 3 分時計と 4 分時計を比較した. 砂粒の動きに注目すると, 3 つの間に特に大きな相異点は見られなかった. また, 砂の総量は表 3 より, １分時計は 3 分時計の約 1/3, 4 分時計の約 1/4 であることが分かる. これより, 砂の総量と砂が落下しきる時間はほぼ比例していることが分かる.
表3

時間の違いによる砂粒の総量の比較

時間

１分

3分

4分

砂の総量

4539.11114094.83
3.1

18648.83
4.1

比率

4.3 考察4.2 節での比較の結果では, 以下のような挙動が見られた.
1. 形状の違いによる落下時間の変化は見られない
2. 傾斜角がゆるやかだと砂粒の流出速度が徐々に低下

3. 砂粒の総量と砂粒が落ちきる時間は比例
1, 3 は, 砂時計の動きとして実際に見られる挙動であり, これらについては, 本研究で作成した砂時計のモデル

本研究では, 粉体の動きを観察できる例として砂時計のシミュレーションをセルオートマトンを用いて行った. 砂の量をセルの高さという情報で扱い, 砂柱の崩壊やそれに伴う蓄積を高さの変化で表した. そして, 落下中の砂の挙動を考慮した上で砂の崩壊, 落下, 蓄積を表現した. セルオートマトンを用いたことで, 容器の形状を変更したときの砂粒の挙動の違いがないこと, 砂粒の量のみによる時間の変化など, 砂時計の挙動を再現できた.
しかし, 容器の側壁の傾斜角の変化によって, 砂粒の崩壊の早さや砂粒が落ちきる時間に変化が発生してしまった. これは, 実際の砂時計では容器の形状によらず砂の量のみで時間が決定されると考えられるので, 実際の砂時計の動きを完全に再現できなかったといえる.
改善すべき点としては, 砂粒が落下する際の砂時計の上部での動きである. 本研究では現在, 表層面での変化のみで砂の崩壊を表現しているが, 本来なら砂粒が落下した場合, その周りの砂粒が入り込むというように, 下部から流れていく. これを解決することにより, 容器の傾斜角によって流出速度が変化することも防ぐことができると考えられる.
また, この砂時計を三次元化することも重要である. 三次元化することにより現実の砂時計に近い挙動を表現できるが, そのためには今まで使用していた一次元オートマトンに代わって, 二次元オートマトンで砂の崩壊, 蓄積のパターンを考慮する必要がある.
さらに, 有限差分法, 有限要素法, SPH 法などの他の手法でのシミュレーションを行い, その結果とセルオートマトンを用いた今回の研究結果との比較・検討をすることも重要な課題である.

参考文献[1] C. P. C. Prado, Z. Olami：Physical Review A, 45,665-669(1992).
[2] 加藤 恭義, 光成 友孝, 築山 洋：“セルオートマトン法―複雑系の自己組織化と超並列処理”, 森山出版株
式会社 1998.
[3] 飯塚 尚彦：“粉体シミュレーション入門”, 産業図書株式会社 1998.
[4] 田中 善弘：“砂時計の七不思議”, 中公新書 1995.
[5] 山 内 千 里：“EGGX / ProCALL version 0.78”,
http://phe.phyas.aichi-edu.ac.jp/~cyamauch/
eggx_procall/.