仮想液体操作―密度変化モデルによる流れの表現
舟橋健司 ∗1
宇佐美健一 ∗2
岩堀祐之 ∗1∗2∗1
名古屋工業大学 情報メディア教育センター
株式会社 メイテツコム
E-mail: kenji@center.nitech.ac.jp

Expression of Liquid Flow Using Density Change Modelin Virtual Liquid Manipulation
Kenji Funahashi∗1 ,Ken-ichi Usami∗2andYuji Iwahori∗1∗1∗2
Nagoya Institute of Technology
and
MEITETSU COM Co.,Ltd.
E-mail: kenji@center.nitech.ac.jp

1

はじめにうな動作を可能としている.
・落下して来る液体を容器で受け止める

バーチャルリアリティ(VR) 技術に対する要望のひとつとして仮想空間中における対話操作, というものが・容器内の液体を別の容器ですくいとる・容器から液体を落下させるある. 従来では, 主に仮想空間に存在する固体 (剛体や本研究では, このモデルの拡張として, 流れの状態に弾性体) を対話的に移動, 変形する研究が行われて来た

ある仮想液体の表現手法の提案を行う [10]. 本手法では[1,2]. しかし, 医療やアミューズメントの分野などにおいて固体だけでなく, 液体や気体に対しても仮想空間で対話操作を行いたいという要望もある.

流れている状態の仮想液体を粒子に基づいて表現する.
粒子個々の挙動を追うのではなく, 一定領域内にある液液体や気体の流れを扱った学問としては流体力学が値を隣接する格子間で変化させる事により液体の移動挙げられる. これは液体などの挙動を厳密に計算するを表現する. なお, このモデルは対話操作を第一の目的には適しているが, 計算量が多いため, リアルタイム処としたものであり, CG によるリアリティの高い映像の理を行う必要がある対話操作には直接適用することは生成や, 挙動の厳密な再現を主な目的とはしていない.

体粒子の密度と平均速度を定義する. そしてこの 2 つの出来ない. また, コンピュータグラフィックス (CG) や

VR の分野においても液体, 気体に関する研究は行われている [3-5]. これらの研究も, リアリティの高い CG を2粒子体積に基づくモデル作成することを目的としており対話操作には向いていない. Stam は流体アニメーション作成の高速な手法を

本節では, これまでに提案した任意形状の容器によ提案している [6]. Navier-Stokes 方程式をもとに, テクる液体の対話操作を実現するための仮想容器, および,スチャマッピングを利用しているこの手法は, 対話的に仮想液体のモデルについて述べる [7-9].

アニメーションを作成可能である. しかし, 入れ物を用

このモデルでは液体を以下の 2 つの状態で考える.

いた液体に対する, すくう, こぼす等の操作への直接的●自由落下 現実では, 空から降って来る雨, 水道の蛇な適用は難しい.

口から流れ落ちる水などがある. ここでは, 液体を粒子我々はこれまでに, 仮想空間における任意な形状の容に基づき表現する. 本研究では微妙な動きよりも対話操器による液体の対話操作を実現する新たな仮想液体操作を目的としているため, 各粒子は大きさを考えない.

作モデルを提案している [7-9]. このモデルでは仮想空

そして, 粒子同士の衝突は起こらず,2 つの粒子が衝突間内に液体との干渉を判定すべき仮想の容器を定義し,して 1 つになるという現象はおこらないものとする. こ

容器と液体との相互作用についてのみ考え, 以下のよれらの粒子は重力と慣性によって動く.

density
43133233.2
velocity

32

容器のモデル

以下では容器内における液体の流れを考える. 液体の流れを考える容器 (以下, 単に容器と呼ぶ) を予め一辺の長さが d , 容積 d3 の格子状に区切っておく. そして区切られた格子それぞれについて, 液体の平均速度を v, 格子内にある液体粒子の密度を D [粒子数/容積]

図 1: 流れている液体のモデルと定義する. この速度と密度を隣接する格子間で変化させるという処理を繰り返し行うことにより, 全体的に●静止状態

現実では, コップの中の水などがある. ここでは, 容器内の対流などを考えず, 完全に静止した液体として扱い, 体積により表現する. すなわち, ある時
刻 t において液体が仮想空間内の容器に存在している
場合, 上述の粒子の性質は一切考慮せず, 容器の中に体積 V (t) の液体が存在していると表現する.
上記の 2 状態間の換算には換算比 N [粒子数/体積] を見れば液体が流れているような動作を実現する. 多くの液体は気体のように圧縮されることはないため, 液面付近以外では液体粒子の密度は一定になるはずである (このときの液体の密度を N と定義する, 格子内の液体の粒子数は N d3 である). しかし, 本手法では液体を隣接する格子間で液体を漸次移動させる方法であるため, 一時的に液体粒子の密度 D が一定の範囲内で増減することを認める.

用いる. 仮に容器内に体積 V の静止状態の水があるとして, その全体積をこぼすと N V 個の粒子と体積 0 の静止状態の水となる. このモデルに基づき作成した実験システムでは, 落下してくる液体を操作者が自由に動かすことのできる仮想の容器で受け止め, 容器の中に溜めることが可能である. また, 操作者が容器を傾けることにより, 容器内の液体をこぼすことも可能である. さらに, 他の容器内の液体をすくう操作も可能である.

3

密度変化モデル

3.3

処理の流れ

ある時刻 t における状態から t + ∆t の状態を求めるには以下に示す処理を行う. その処理を繰り返し行う
ことにより液体の流れを表現する. 以降, それぞれの処理について説明する.

(1) 密度の変化の実行(1.a) 各格子内の密度の変化量の計算 まず始めに, 時
刻 t から t + ∆t に状態が変化する際に各格子から他の格子へ移動する液体の粒子数を計算する. 考え方は,格子内に一様分布している液体粒子を平均速度分だけ3.1

液体の流れの表現方法

移動させる, というものである. 時刻 t において, 座標(i, j, k) の位置にある格子 L(i,j,k) 内の液体粒子の密度流れている状態の液体は現実では, 川の水や雨樋をを D(t)(i,j,k), 平均速度を v(t)(i,j,k) と定義した時, 移動流れる水などがある. 液体を粒子に基づく方法で表現する液体の密度は以下の処理を行うことにより求まる.

しようと考えた場合 (図 1 左), 液体粒子それぞれの位置や速度をもとに全体の液体の挙動を考える必要があ

step1 立方体の想定 格子 L(i,j,k) 内の液体が密度 N でD(t)(i,j,k) d3あったときの体積はであらわせる. 次に, 格Nり, 計算量が多くなる. そのため, リアルタイム処理を子 L(i,j,k) と重心が同じ座標にあり, 体積が考える際には個々の液体粒子の挙動を考えるという手の立方体を定義する (図 2(a) 点線部分). この立方体は法を直接用いることはできない.

格子内にある液体を想定するものである. これを次のス

D(t)(i,j,k) d3
N

そこで, 前述した 2 つの状態の中間的な状態を考え

テップ以降で変形, 移動した後の位置を基準に L(i,j,k)る. まず液体の流れを考慮する部分を予め格子状に区から隣接する格子に移動する液体の量を決定していく.

切っておく. この格子内の液体粒子それぞれの挙動をなお, 簡単のために 2 次元の図で表現しているが, 実際見るかわりに, この領域内では液体粒子が一様に分布しの処理は 3 次元で行う.ていると仮定して, 一定領域内に含まれる液体粒子の個の 2 つの値を隣接する格子間で変化させることによっstep2 圧力を考慮した形状の変形 注目している格子内の液体密度に対して隣接格子内の液体密度が高い場
合, 圧力により注目格子の外へ移動する液体は少なくて液体の流れを表現する.

なる. 逆に注目格子に対して隣接格子内の液体密度が数, つまり液体の密度と平均速度を考える (図 1 右). こ低い場合には外へ移動する液体は多くなる. このこと

D(t)(i,j,k) d3
NYを考慮し step1 で求めた立方体の形状を変形する. 立(i,j,k)

方体の各面の変化量は D(t)Y(−1,0,0)(0,1,0)m(i,j,k)m(i,j,k)(1,0,0)m(i,j,k)と, 周囲の格子内の(0,−1,0)m(i,j,k)

液体の密度 D(t)(i+p,j+q,k+r) (ただし,step2 においては (p, q, r) = (−1, 0, 0), (1, 0, 0), (0, −1, 0), (0, 1, 0),(0, 0, −1), (0, 0, 1) の, 面で隣接する 6 通り) をそれぞ(p,q,r)
れ比較し, その差に比例する距離 m(i,j,k) だけ立方体の各面の位置を変化させる (図 2(b) 破線部分).

(a)

X(i,j,k)(i−1,j −1,k)∆tv(t)(i,j,k)

YX(b)B(i,j,k)Y B(i,j,k)(p,q,r)m(i,j,k) = ρ(D(t)(i,j,k) − D(t)(i+p,j+q,k+r) )ρ : 比例定数step3 直方体の移動 次に step2 によって変形された直方体を格子内の液体の平均速度 v(t)(i,j,k) を用いて移動させる (図 2(c) 太線部分). これは時刻 t において(c)格子 L(i,j,k) にあった液体が流れる様子を表している.
直方体の移動距離は ∆tv(t)(i,j,k) で表せる.

step4 移動する密度の計算

XX(d)図 2: 格子内密度の変更量の計算

step3 の処理により移動する液体 (粒子) の割合となる. 直方体のうち,L(i,j,k)(2) 速度の変更(2.a) 密度の変化に伴う速度の変更 密度の更新を行ったことに伴い, 各格子内の液体の平均速度も変更するの近傍の格子内部に含まれる部分の体積を, それぞれ必要がある. ある格子 L(i,j,k) において, 密度の更新をした直方体が, どの格子の中にどれだけの割合で含ま(i,j,k)

れているかを計算する. この割合が, Lから移動(i+p,j+q,k+r)

B(i,j,k)とする (図 2(d) 斜線部分). ただし, step4においては p, q, r のとる値はそれぞれ −1, 0, 1 のいずれかである. 直方体のうち, 格子 L(i,j,k) の内部に含ま(i,j,k)

れている部分の体積は B(i,j,k) であらわせる. このようにおいた時, L(i,j,k) から単位体積の格子 L(i+p,j+q,k+r)(i+p,j+q,k+r)に移動する液体の粒子数 T(i,j,k)は以下の式で求まる. 他の格子へ移動せず元の格子に留まる液体は(i,j,k)T(i,j,k) , すなわち p = q = r = 0 により求まる.
(i+p,j+q,k+r)(i+p,j+q,k+r)

T(i,j,k)= D(t)(i,j,k) 1

B(i,j,k)(i+u,j+v,k+w)
u,v,w=−1 B(i,j,k)

ただし, 液体の移動先の格子内の密度が設定されたしきい値よりも高い場合は, その格子へ液体が移動して更に(i+p,j+q,k+r)=0と密度が高くなるのを防ぐために, B(i,j,k)

して計算する. また, 直方体の移動先が容器の壁面より(i+p,j+q,k+r)も外部に出てしまった場合も同様に B(i,j,k)=0として計算する.
以上の step1 から step4 までの処理を全ての格子について行うことにより, 液体の移動量を計算する.

(1.b) 全格子の液体密度の変更 前述の方法により各格子について密度の変化量を計算した後, 全ての格子の密度を更新する. ある格子 L(i,j,k) に注目した時, 液体を移動させた後の密度は次式で示される.

D(t + ∆t)(i,j,k) =1
p,q,r=−1(i,j,k)T(i+p,j+q,k+r)(i,j,k)行った後の平均速度 v(t)Tは以下の通りである.

1vT =(i+u,j+v,k+w) (i,j,k)T(i+u,j+v,k+w)
u,v,w=−1 v(t)
1(i,j,k)p,q,r=−1 T(i+p,j+q,k+r)(2.b) 液体の加速 液体の流れを考えている容器が移動する際, 容器内の液体は元の位置に留まろうとするので, 容器に加えられた力とは逆向きに反作用の力がかかる. 時刻 t における容器の中央の座標を C(t) とすると, 液体にかかる加速度 a(t) はa(t) = −C(t) − 2C(t − ∆t) + C(t − 2∆t)∆t2である. また, 液体には常に重力加速度 G がかかっている. これらにより, 格子 L(i,j,k) 内の液体の速度は以下のようになる.

v(t)(i,j,k)= v(t)T + ∆tG + ∆ta(t)a(3) 例外処理(3.a) 壁面における液体の反射 容器の壁面 (底面) 部分においては, 液体の速度の壁面に垂直な方向成分は反対方向へ反射する. 本手法では液体を減衰率 η で非弾性衝突させる (図 3(a)). ただし振動を抑制するために,壁面に対して垂直な方向の速度成分が一定のしきい値

f 以下の時は反射は行わず, 垂直方向の成分を 0 とする(i,j,k)(図 3(b)). 格子 L(i,j,k) の速度が v(t)aの時, 壁面

(i,j,k)v(t)h(i,j,k)v(t)rた後の速度 v(t + ∆t)(i,j,k) は次式で表せる.

(i,j,k)v(t)r(i,j,k)−ηv(t)p(i,j,k)v(t)p(i,j,k)v(t)a(i,j,k)v(t)p(i,j,k)v(t)a(i,j,k)v(t)r|vmax |
(i,j,k)|v (t)r|(i,j,k)(|v(t)r| > |vmax |)(i,j,k)(|v(t)r| ≤ |vmax |)容器の大きさを考慮して格子の大きさを設定することにより, 一般に扱う液体に関しては問題は小さいと考(b)(a)
図 3: 壁面部分における液体の反射

えられる.

3.4
(i,j,k)に対して垂直な方向の速度成分を v(t)p, 平行な方(i,j,k)とすると, 壁面での反射を考慮し向の成分を v(t)h(i,j,k)た後の速度 v(t)rは以下のようになる.

(i,j,k)(i,j,k)− ηv(t)p v(t)h(i,j,k)vr =v(t)h
(i,j,k)v(t)a(i,j,k)v(t)rv=(i,j,k)(f < v(t)p)(i,j,k)(f ≥ v(t)p)(壁面以外)(3.b) 容器の底面付近での処理 容器の底面付近の格自由落下状態の液体と容器

ここでは自由落下状態の液体が容器に入り, 流れの状態になる時, 逆に流れの状態から自由落下状態の液体になる時の遷移モデルについて説明する.
自由落下状態から流れの状態へ

自由落下状態の液体(ある粒子) が容器の内部にあった場合, その液体を粒子としては消滅させ, その粒子の存在した座標を含む格子内の液体の密度を増加させる. また速度の更新時に粒子の最終速度を考慮する.
流れの状態から自由落下状態へ

容器の口を構成する子では, 重力により下向きの速度をもった液体が常に流格子内の液体の移動先が, 容器の外部であった場合, 移れ込んでくる. そのため, 密度が高い状態が続いてしま動した液体の密度に比例した量だけ自由落下状態の粒い, 時間が経過しても容器内の液体の密度は均等には成子を作成する.

らないという状況が生じる. この問題を改善するため,容器の底面の方では例外的な処理を行う.
条件 1: 格子内の密度がしきい値より大きい.
条件 2:

格子の下方 3 方向が容器の壁面であるか, ここで示す例外処理が行われている状態であり, 下方へと液体が流れてはいけない状況である.

4

実験および結果
上述のモデルに基づいた仮想液体操作のための実験システムを, グラフィックワークステーション SGI OC-以上の条件を満足した格子について, 格子の各面を容TANE2 上で C 言語により実現した. 液体の流れを考える仮想容器は操作者が自由に移動することが可能で,容器の位置と傾きは低周波の磁界を利用したモーショ器の壁面と同様の扱いにする. つまり, 隣接する他の格ンセンサにより入力される. 液体の流れを考慮する容子からこの格子へ移動してこようとしてくる液体につ器の形状は縦横 20cm, 高さ 6cm を想定している. まいては, 移動させず, 平均速度についてはこの格子の面た, 格子の一辺の長さは 1cm, 容積 1cc とする. 自由落に垂直な方向に反射させる. ただし, この格子から他の下状態の液体粒子は 1 個 0.1cc として扱う. 限られた格子へ移動する液体については通常どおり移動させる.

画面内で実験を行うために, 重力加速度を 1m/s2 に設このような処理を条件 1,2 を満たさなくなるまでの間,定する.

条件 3: 1,2 を満たす状況が一定時間継続している.

つまり格子内の密度が一定値以下になるか, 下方へと液本システムを使用して実際に操作している例を 図 4体が流れて行ける状態になるまで続ける.に示す. 体験者がセンサの取り付けられたカップを操作(3.c) 速度の制限 最後に, 液体の速度が限界値を超えすることにより, 仮想容器で仮想液体を操作可能なことないように調整を行う. このモデルでは速度 v は一定がわかる. 操作時の画面全体は 図 5 に示す通りである.

速度 |vmax | 以上にはならないものとしている. これはこの時, 描画更新速度は毎秒 20 frame 程度であった.

v が大きくなりすぎると, 単位時間の間に隣接する格子よりも遠方まで液体が移動することを考慮する必要があるためである. つまり, 液体の最高速度 |vmax | は時
間 ∆t の間に格子の 1 辺の長さ d だけ進む速度とする.
(i,j,k)

格子 L(i,j,k)の速度が v(t)rの時, 調整をおこなっ

次に, 液体の流れを考えている容器をすこし左側に傾けた時の様子を図 6 に示す. 容器を傾けたことによって液体が左側へと流れて行き, 容器から外に出た液体が自由落下状態の粒子となっている. 液体が少ない時の様子を 図 7 に示す. ここでは容器が左手前方向に傾いていて, 少量の液体が溜まっている状態から容器を右

参考文献

奥方向に傾けた時の液体の動作を示している. 左前方に溜まっていた液体が広がりながら流れて行き, 最終的に右奥方向に向かっていく様子が分かる. 実際に研究室内の数名の学生に体験実験をしてもらったところ,本当に液体を扱っているように感じる, 液体が流れているように見える, などといった積極的な意見を頂いた.

[1] 舟橋健司, 安田孝美, 横井茂樹, 鳥脇純一郎: “仮想空間における仮想手による道具操作のための知識とモデル”, 日本バーチャルリアリティ学会論文誌,Vol.3, No.3, pp.167-176 (1998)

このモデルの現状の問題点としては, 液体の量が少[2] 藤本英雄, 佐野明人, 磯部玄, 胡摩心一郎, 竹内博
良: “マイクロドームシステムを用いた力覚縫合ない時, 液面が完全に水平にはならない, ということがシミュレータ”, SICE SI2001 講演論文集, pp.405-挙げられる. 例として, 容器の底面が水平方向と平行で,406 (2001)

容器の底面にある格子内の液体の密度が低く, 各格子内の液体の速度が全て 0 であるような状況を仮定する.
この時, 密度変化を考える際に使用する直方体 (図 2(c)の太線部分にあたる直方体) がもとの格子の内部に含まれてしまい, この格子内の液体は他の格子へ移らない. そのため, 液面が水平にならない状態のままで安定してしまうという問題が存在する.

5むすび
仮想空間における流れの状態にある液体の表現を, 液体の流れを考える容器を格子状に区切り, 局所的に一様に分布する粒子を想定した密度変化モデルを用いて[3] Michael Kass, Gavin Miller:“Rapid, Stable

Fluid Dynamics for Computer Graphics,” Proc.
SIGGRAPH 90, pp.49-57 (1990)[4] N. Chiba,

S. Sanakanishi,

K. Yokoyama,

I. Ootawara: “Visual Simulation of Water Currents Using a Paticlebased Behavioural Model,”The Journal of Visualization and Computer Animation, Vol.6, pp155-171(1995)[5] Jim X. Chen, Niels da Vitoria Lobo, Charles E.
Hughes, J. Michael Moshell: “Real-Time Fluid
Simulation in a Dynamic Virtual Environment,”
IEEE CG&A, Vol.17, No.3, pp.52-61 (1997)実現した. 本システムでは容器を傾けることによって流れている状態の液体を自由に操作することができ, さらに自由落下してくる液体を受け止めること, 逆に流れている状態から液体をこぼして自由落下状態にすることも可能である.
今後の課題として, 前述の液面が完全に水平にはならない問題に対応する必要がある. 本システムでは経[6] Jos Stam: “Stable Fluids,” Proc. SIGGRAPH
99, pp.121-128 (1999)[7] 舟橋 健司, 岩堀 祐之: “仮想容器による仮想液体の対話操作モデルと一実現法”, 日本バーチャルリ
アリティ学会論文誌, Vol.5, No.4, pp.1087-1094(2000)

験的に設定した格子の大きさや液体の粒子との換算比に関しても, 実験, 考察しなければならない. また実際

[8] Kenji Funahashi, Yuji Iwahori: “Virtual Liq-の液体では, 非圧縮性流体の場合には, 圧力は密度, 重uid Manipulation Using General Shape Vessel,”
Proc. IEEE VR2001, pp.283-284 (2001)

力と深さの積により考えることができる. 運動と位置,圧力により伝えられるエネルギーを考慮し, 物理法則

[9] Kenji Funahashi, Yuji Iwahori: “Representationに基づくパラメータを導入した新たなモデルへと発展

of Swinging Liquid on Virtual Liquid Manipulation,” Proc. ICAT2001, pp.179-184 (2001)

させることにより, さらに自然な流れを表現する必要がある. 液体の挙動だけでなく, 臨場感を高めるためにCG による表現方法の向上や波の音などの表現についても取り組みたい.

[10] 宇佐美健一, 舟橋健司, 岩堀祐之, 岩田彰: “格子モデルを用いた流れる液体に対する対話操作の検討”, SICE SI2001 講演論文集, pp.403-404 (2001)

謝辞
最後に, 本研究を進めるにあたり有益な議論を頂いた本研究室諸氏, および卒業生に感謝する.

図 4: 実際の操作の様子

図 5: 操作時の画面表示

図 6: 実験結果の連続画像

図 7: 液体が少ないときの動作