2008 年度 修士学位論文
トポロジ変化の実時間レンダリング
2009 年 2 月 13 日
大阪大学 大学院基礎工学研究科
機能創成専攻 生体工学領域
荒井 良祐
主査: 日付: 副査: 日付: 概要
物体は,原子と原子間の結合で定義されるネットワークモデルであり,接続形態であるトポロ
ジを有している.物体に大きな外力が作用し,生じた応力が限界値を超えると破壊が生じる.破
壊により原子間の結合が切れ,物体のトポロジが変化する.現実世界では,破壊操作を繰り返し
行うことが困難であるが,計算機上の仮想空間では,仮に破壊現象を扱うモデリング手法が存在
すれば,破壊操作が可能となる.近年,様々な分野で力覚提示可能な VR システムが開発され,医
療や芸術等の分野で用いられている.変形物体を扱う多くのシステムでは,有限要素モデルが用
いられているが,トポロジ変化を伴う破壊現象を可能とするモデリング手法は十分に確立してい
ない.主要因は,トポロジ変化時に必要となる幾何形状,および,剛性マトリクスの再構築計算
コストが高いことである.トポロジ変化時に,物体全体ではなく局所部分の再構築計算のみを行
うことで,計算時間を削減できると考えられる.
本研究では,ユーザ操作時に仮想物体の応力状態に基づいた破壊現象を実時間処理し,力覚提
示可能とする手法を提案する.提案手法は,破壊によるトポロジ変化時の計算コストを削減する
手法,ならびに,トポロジ変化を詳細に反映した反力を算出する手法で構成される.計算コスト
削減は,き裂先端領域近傍の局所領域から構成されるローカルオブジェクトを定義し,ローカル
オブジェクトについてのみ再構築計算を行うことで実現する.破壊進行に伴い,ローカルオブジェ
クトは,き裂先端近傍を囲むように更新される.また,反力算出のために,前処理段階で多様な
トポロジとユーザ操作における反力を記録し,実時間処理中には記録した反力を補正する.さら
に,前処理で用いるパラメータは,ヒトの力覚的分解能を満たす程度に限定できると考え,被験
者実験によってヒトの破壊進行角度に対する力覚的分解能を調査し,パラメータを限定した.
提案手法により,仮想物体の応力状態に基づいた破壊現象の実時間処理,および,力覚応答が
可能であることが確認された.ユーザは仮想物体を引きちぎる操作が可能となった.また,従来の
ように物体全体の再構築計算を行った場合に比べて,提案手法では計算時間を削減することがで
きた.さらに,仮想物体の節点数を変えて再構築時間を比較した結果,節点数の増加に伴う計算
時間の増加が,従来法と比較して小さかった.要因は,提案手法では,従来法と比べて再構築計
算の対象となる節点数が少ないこと,また,再構築計算時間時に,節点数の三乗のオーダの計算
コストを要する逆行列計算を行わなかったことが挙げられる.さらに,反力補正の有効性を示す
ために,補正を行う場合と補正を行わない場合の誤差を比較した.結果,比較した多くのパター
ンにおいて,補正を行った場合の方が,補正を行わない場合に比べて誤差が小さくなった.また,
誤差の程度についても,補正を行った場合の方が,補正を行わない場合に比べて誤差の平均値,標
準偏差ともに小さくなった.以上の結果から,提案した手法によって,再構築計算時間の削減,お
よび,反力誤差の減少が確認され,提案手法は有効であると言える.
キーワード:力覚提示,破壊現象,トポロジ変化,有限要素法,実時間計算
Abstract
An object can be defined by a network model consisting of atoms and connections among atoms.
An object has a topology as a pattern of connections. When stress caused by large external
forces exceeds the limit, fracture phenomena occurs. Connections among atoms are broken and
a topology changes by fracture. In virtual world, if the modeling method of treating fracture
phenomena exists, an object can be repeatedly fractured without spoiling the object physically
in real world. In the past decade, VR systems that can present reaction forces are intensively
studied and applied to the various fields such as medical care and art. A lot of recently developed systems with deformable objects use the finite element models. However, the technique for
topological change by fracture phenomena in real-time has not been established. One of the key
factors is that the calculation cost of restructuring of the mesh and inverse stiffness matrix is
high. The calculation time of restructuring could be reduced by restructuring only local region
surrounding a fracture region instead of global object.
This study proposes methods to treat a fracture phenomena based on stress condition with
user’s interactive manipulations in real-time. Proposed method consists of two techniques. One
of the techniques reduces restructuring calculation costs in topological changes by fracture and
the other techniques present reaction forces that reflects topological change of objects to the
user. In order to reduce calculation cost of restructuring, local object is defined as a small
region, in only which restructuring calculation is applied, surrounding a crack. As a crack propagates, local object is updated surrounding a new crack. In order to calculate reaction forces,
this study proposes a recording and replaying methods with interpolation of crack orientations.
Moreover, in order to limit the parameters of pre-calculated topological patterns to extent in
which haptic resolution of human is fulfilled, the parameters are determined with an experiment
that examines haptic resolution of human for the orientation of fracture progress.
The developed system showed that the proposed model enabled to treat a fracture phenomena
based on stress condition with user’s interactive manipulations in real-time. The user can pull
and tear a virtual object in real-time. The calculation time of restructuring with the proposed
method became shorter than that one by conventional method. Moreover, the difference of
calculation time between conventional and proposed method has grown as the number of nodes
increases. One of the factors is that the number of calculated nodes in the proposed methods
was smaller than those of the conventional method. The other factor is that inverse matrix
calculation which requires the calculation cost of cube order of model size was unnecessary.
In order to examine the effect of proposed interpolation method of reaction force, errors with
interpolation were compared with the ones without interpolation. As a result, the errors with
interpolation were smaller than the errors without interpolation in many patterns. Moreover
means and standard deviation with interpolation were smaller than the ones without interpolation. From the above-mentioned results, a decrease of calculation time of restructuring and a
decrease of error were confirmed by the proposed method, so the effect of proposed method was
shown.
Keyword : Haptic representation, Fracture phenomenon, Topology change, Finite Element
Method, Real-time computation
目次
第 1 章 はじめに
1
第2章
2.1
2.2
2.3
固体の破壊現象
固体の材料特性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
弾性体 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
破壊現象 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
5
12
第3章
3.1
3.2
3.3
力覚提示技術と応用
有限要素モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
力覚提示デバイス . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
アプリケーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
18
25
29
第4章
4.1
4.2
4.3
実時間トポロジ変化モデル
提案手法概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
再構築計算時間削減 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
トポロジ変化を詳細に反映した反力提示 . . . . . . . . . . . . . . . . . . . . . . .
33
33
36
42
第5章
5.1
5.2
5.3
5.4
システム実装
システム概要 . .
提案モデル実装 .
破壊操作の結果例
定量評価 . . . . .
47
47
48
54
55
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
第 6 章 おわりに
58
謝辞
59
業績
60
参考文献
61
図目次
1.1
破壊現象に伴うトポロジ変化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12
2.13
2.14
外力と変位の関係 . . . . . . . . . . . . . . . .
本研究で仮定する外力と変位の関係 . . . . . .
応力の定義 . . . . . . . . . . . . . . . . . . .
力のモーメントの釣り合い . . . . . . . . . . .
主応力と座標変換 . . . . . . . . . . . . . . . .
変形状態 . . . . . . . . . . . . . . . . . . . . .
xy 面内の変形状態 . . . . . . . . . . . . . . .
破壊の分類 . . . . . . . . . . . . . . . . . . .
最大主応力条件による破壊の概念図 . . . . . .
最大せん断応力条件による破壊の概念図 . . .
せん断ひずみエネルギ条件による破壊の概念図
長さ 2a のき裂を有する弾性板 . . . . . . . . .
代表的な三つの負荷条件 . . . . . . . . . . . .
き裂先端を原点とした極座標系 . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3
4
5
6
8
9
9
12
13
14
14
15
15
16
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17
有限要素法の概念図 . . . . . . . . . . . . . . .
固定節点,自由節点,接触節点の概念図 . . . .
面上の任意の点に接触した場合の反力計算 . . .
面積座標の例 . . . . . . . . . . . . . . . . . . .
PHANToM Omni の概観 [43] . . . . . . . . . .
PHANToM Premium の概観 [43] . . . . . . . .
Falcon の概観 [44] . . . . . . . . . . . . . . . . .
SPIDAR の概観 [45] . . . . . . . . . . . . . . .
HapticMaster の概観 [46] . . . . . . . . . . . . .
CyberGrasp と CyberForce の概観 [47] . . . . .
Gravity Grabber[48] . . . . . . . . . . . . . . .
流体グローブ [49] . . . . . . . . . . . . . . . . .
仮想切断シミュレータ [50] . . . . . . . . . . . .
剥離シミュレータ [55] . . . . . . . . . . . . . .
Haptic texture の概念図 [56] . . . . . . . . . . .
記録・再生型アプローチを用いた変形モデル [58]
Haptic video の概念図 [59] . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
19
22
24
24
25
25
26
26
27
27
28
28
29
30
30
31
32
4.1
4.2
提案手法で用いるモデルの概略図 . . . . . . . . . . . . . . . . . . . . . . . . . . .
処理の流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
34
1
1
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
4.11
入力を受けて応力が算出された状態 . . . . . . . . . . . . . . .
破壊進行領域決定 . . . . . . . . . . . . . . . . . . . . . . . . .
破壊進行 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ローカルオブジェクト更新 . . . . . . . . . . . . . . . . . . . .
交差する面と線分 . . . . . . . . . . . . . . . . . . . . . . . . .
交差線分算出 . . . . . . . . . . . . . . . . . . . . . . . . . . .
新たなき裂形状形成の様子 . . . . . . . . . . . . . . . . . . . .
一回目の現在の破壊進行角度と前処理した破壊進行角度の関係
二回目の現在の破壊進行角度と前処理した破壊進行角度の関係
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
37
38
38
39
40
41
41
43
45
5.1
5.2
5.3
5.4
5.5
5.6
5.7
システム全体図 . . . . . . . . . . . . .
モデル外観 . . . . . . . . . . . . . . .
被験者実験の様子 . . . . . . . . . . . .
被験者実験に用いた見本とモデル A,B
クラスタ外観 . . . . . . . . . . . . . .
ユーザ操作による仮想物体の破壊結果 .
計算時間測定結果 . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
47
48
50
51
53
54
55
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
表目次
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
機器の仕様 . . . . . . . . . . . . . . . . . . .
実装環境 . . . . . . . . . . . . . . . . . . . . .
実装に用いたモデルの仕様 . . . . . . . . . . .
被験者実験の結果 . . . . . . . . . . . . . . . .
前処理パラメータ . . . . . . . . . . . . . . . .
クラスタ計算機の仕様 . . . . . . . . . . . . .
再構築時間比較に用いたモデルの仕様 . . . . .
提示反力誤差を調査する際に用いたパラメータ
提示反力評価結果(進展回数 1 回) . . . . . .
提示反力評価結果(進展回数 2 回) . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
48
48
49
52
52
53
55
56
56
56
第 1 章 はじめに
私たちは日々の生活の中で,割り箸を割る,もしくは,お菓子の袋を開けるというように,身の
回りの物に対して破壊現象 (Fracture phenomenon) を含む操作を行っている.普段意識すること
は少ないが,破壊とは一体どのような現象であるかを考えてみたいと思う.多くの物体は,外力が
作用すると変位が生じて変形するが,変位が微小であれば弾性変形をする.弾性とは,操作によ
り外力を作用しているときは変形し,操作を止めると元の状態に戻る性質である.次第に外力を
大きくしていくと,局所的に応力が集中し,やがて破壊が生じる.微視的には,物体を構成してい
る原子 (Atom) 間の結合 (Connection) が切れ,分離する現象である.一方で,物体を原子と原子
間の結合で構成されるネットワークと見ると,物体は原子の接続形態であるトポロジ (Topology)
を有していると言える.破壊により原子間の結合が切れることは,ネットワークの構成要素であ
る結合が切れて,トポロジが変化することと同義である [1].つまり,図 1.1 に示すように,破壊
現象によって物体のトポロジが変化する.
Connections Atoms
Object
Microscopically vision
Fracture (Topological change)
Microscopically vision
Connections are broken
図 1.1: 破壊現象に伴うトポロジ変化
破壊現象が不可逆な性質を持つために,破壊を伴う操作を繰り返し同じ状況で行うことは,現
実世界では困難である.しかし,計算機によって作られた仮想空間内では,仮想物体に対して操
作を行うため,何度でも繰り返し同じ操作を行うことができる.仮に,仮想空間で破壊現象が実
現可能なモデルを構築可能となれば,現実世界では困難であった,破壊現象を繰り返し行うこと
1
ができるのではないかと考えられる.
近年,Virtual Reality(VR) に関する研究が盛んに行われている. VR とは,物体や現象の本質
的な効果を人間の五感に提示することで,実際に物体が存在していない,もしくは現象が起きて
いない状況下において,物体や現象の存在を感じさせる技術である. VR 研究では,物理現象や
物体操作の模擬が行われており,手術手技の訓練など医療分野 [2]-[25],芸術分野 [26]-[28] やエン
ターテイメント [29]-[33] で応用されている.VR 環境下における物理変形モデルとして,有限要
素モデルが数多くの研究で用いられている.有限要素モデルでは,物理計算に有限要素法 (Finite
element method : FEM)[34]-[36] を用いることで,ユーザ操作時における物体の変形や反力を,連
続体力学に基づいて高精度に模擬する事ができる.また,バネ質点モデルや充填球モデルと比較
して,応力とひずみといった物理量が算出できることや,ヤング率やポアソン比といった材料定
数をそのまま用いることができ,現実の物体を容易にモデリングできるといった利点がある.現
状では,弾性変形の実時間模擬が可能となっているが,破壊現象の実時間模擬は困難である.主
要因は,破壊によって変化する剛性マトリクスの逆行列の再構築計算コストが高いことである.
本研究では,従来の有限要素モデルでは困難であった,破壊現象の実時間処理を可能とする手
法を提案する.仮想物体としては,等方性の線形弾性体を用い,初期状態でき裂を有している板
状の三次元モデルを用いる.再構築計算コスト削減は,物体全体ではなく局所領域についてのみ
再構築計算を行うことで実現する.再構築計算を行う局所領域は,破壊力学 [37] の知見から,き
裂先端領域近傍に限定できる.また,トポロジ変化を詳細に反映した反力提示を行うために,前
処理段階で,多様なトポロジと複数の操作における反力データを記録し,実時間処理中には物体
のトポロジとユーザの操作に応じて補正して再生する.
2
第 2 章 固体の破壊現象
固体の破壊現象を扱う上で,固体の力学的な知見である材料特性や他の分野における破壊現象
の取扱い方を学ぶことは,必要不可欠である.本章では,初めに固体の材料特性について述べた
後,弾性体の性質,理論について記す.次に,破壊現象を扱うための基礎知識として,破壊の種
類や破壊基準,および,破壊現象を定量的に取り扱う破壊力学の知見について記す.
2.1
固体の材料特性
Force
物体に徐々に大きな外力を作用させる場合を考える.作用させる外力が小さい間は,物体は弾
性変形し,徐々に外力を大きくしていくと,物体は塑性変形を生じる.塑性とは,外力を取り除
いても完全に元に戻らず,変位が残る性質である.さらに,大きな外力を作用させると,やがて
物体は破断する.図 2.1 に,物体の外力と変位の関係の模式図を示す.ただし,外力を縦軸,変位
を横軸に示す.
Fracture point
Yield point
A : Linear elastic region
B : Non linear elastic region
C : Plastic region
A
B
C
Displacement
図 2.1: 外力と変位の関係
図 2.1 において,弾性変形を示す区間 A,および,区間 B を弾性域 (Elastic region),塑性変
形を示す区間 C を塑性域 (Plastic region) と呼ぶ.また,弾性域と塑性域の境界を降伏点 (Yield
point) と呼ぶ.さらに,弾性域は,変位が微小な場合は,変位と外力が線形関係にあり,線形弾
性域 (Linear elastic region) と呼ぶ.変位が線形弾性域を超え,変位と外力が非線形関係にある区
間を,非線形弾性域 (Non linear elastic region) と呼ぶ.また,破断を示す点を破断点 (Fracture
point) と呼ぶ.
3
Force
物体は,破断に至るまでに,線形弾性域,非線形弾性域,塑性域という段階を経る.しかし,本
研究では,破壊現象によるトポロジ変化時の実時間処理と力覚提示を主目的とするため,簡単の
ために線形弾性体を取り扱う.線形弾性体では,非線形弾性域と塑性域が微小であり,線形弾性
域の直後に破断すると仮定する.図 2.2 に,本研究で仮定する,物体の外力と変位の関係を示す.
Fracture point
A : Linear elastic region
A
Displacement
図 2.2: 本研究で仮定する外力と変位の関係
4
2.2
弾性体
本節では,弾性体の中でも応力とひずみ間の関係が線形である線形弾性体について,基本的な
力学量である応力とひずみ,および,両者の関係について説明する.
2.2.1
応力
弾性体に外力が作用する時,弾性体内部には作用している外力に抵抗する力(内力)が生じる.
内力を単位面積当たりの力として表したものが応力である.応力は面に垂直にはたらく垂直応力
と面に沿ってはたらくせん断応力とに分別される.図 2.3 に,物体内部の微小領域における応力の
定義を示す.σij は i 軸に垂直な面にはたらく j 方向の垂直応力であり,τij は i 軸に垂直にはたら
く j 方向のせん断応力である.
External force
Object
Microscopical region
図 2.3: 応力の定義
5
図 2.3 において,Δx, Δy, Δz は微小量であるとして,式 (2.1) に示す近似を行う.
∂σxx
Δx
∂x
∂σyy
Δy
σyy +
∂y
∂σzz
Δz
σzz +
∂z
∂τxy
τxy +
Δx
∂x
∂τyx
Δy
τyx +
∂y
∂τzx
Δz
τzx +
∂z
∂τxz
Δx
τxz +
∂x
∂τyz
τyz +
Δy
∂y
∂τzy
Δz
τzy +
∂z
σxx ≈ σxx +
σyy ≈
σzz ≈
τxy ≈
τyx ≈
τzx ≈
τxz ≈
τyz ≈
τzy ≈
応力には九つの成分が存在し,応力テンソルとして式 (2.2) のように表される.
⎡
⎤
σxx τxy τxz
⎢
⎥
⎣ τyx σyy τyz ⎦
τzx τzy σzz
(2.1)
(2.2)
また,力のモーメントの釣り合いについて考える.図 2.4 のように z 軸に垂直な面にはたらく応
力によるモーメントを示す.
P axis
図 2.4: 力のモーメントの釣り合い
6
軸 P についての,力のモーメントの釣り合いは式 (2.3) のようになる.
Δy
Δx
Δy
Δx
Δy
+ σyy ΔxΔz
+ τzx ΔxΔy
+ τzy ΔxΔy
− σxx ΔyΔz
2
2
2
2
2
Δx
Δy
Δx
− τzx ΔxΔy
− τzy ΔxΔy
+ τxy ΔyΔzΔx − τyx ΔxΔzΔy = 0
−σyy ΔxΔz
2
2
2
σxx ΔyΔz
(2.3)
式 (2.3) の左辺では,第 1 項目から第 8 項目までが互いに打ち消しあい,結果として式 (2.4) が
成立する.
τxy = τyx
(2.4)
同様に x 軸,y 軸について垂直な面について考えると,式 (2.5),式 (2.6) が成り立つ.
τyz = τzy
(2.5)
τzx = τxz
(2.6)
最終的に,九つの応力成分の中で独立な成分は σxx , σyy , σzz ,τxy , τyz , τzx の六つのみとなる.ま
た,六つの独立成分を用いて応力ベクトル σ は式 (2.7) のように定義される.
⎡
⎤
σxx
⎢
⎥
⎢ σyy ⎥
⎢
⎥
⎢ σzz ⎥
⎥
σ=⎢
(2.7)
⎢ τ ⎥
⎢ xy ⎥
⎢
⎥
⎣ τyz ⎦
τzx
また,応力は式 (2.2) で表されるように二階のテンソル量であり,考える座標系によって値が変
わる.せん断応力が全て零になるように座標系を設定したとき,最大応力成分を主応力と呼び,座
標系の基底ベクトルを主軸と呼ぶ.主応力や主軸は,式 (2.2) の応力テンソルを対角化することで
求めることができる.応力テンソルを対角化したときの三つの固有値が主応力,三つの固有ベク
トルが主軸となる.式 (2.8) に,応力テンソルと主応力,主軸の関係を示す.なお,三つの主応力
(固有値)を σ1 , σ2 , σ3 とする.また,対応する固有ベクトルを v 1 , v 2 , v 3 とし,固有ベクトルから
構成される行列を V = [v 1 , v 2 , v 3 ] とする.
⎡
⎤
⎡
⎤
σxx τxy τxz
σ1 0 0
⎢
⎥
⎢
⎥
V −1 ⎣ τyx σyy τyz ⎦ V = ⎣ 0 σ2 0 ⎦
τzx τzy σzz
0 0 σ3
(2.8)
最大となる主応力を最大主応力,最小となる主応力を最小主応力,残りの一つを中間主応力と
いう.また,固有ベクトルから構成される行列 V が,三つの主軸で構成される座標系への変換行
列となる.図 2.5 に,主応力と座標変換の概念図を示す.
7
Before transformation
After transformation
All shear stress are zero
Object
Resultant force
図 2.5: 主応力と座標変換
8
2.2.2
ひずみ
弾性体内に応力が発生すると弾性体は変形する.変形の度合を示すために,ひずみと呼ばれる
変位の無次元化量を導入する.物体の点 P が変形によって点 P に移動する変形状態を考える.図
2.6 に,想定する変形状態を示す.
図 2.6: 変形状態
点 P から x 軸に沿って Δx だけ離れた点を点 X, y 軸に沿って Δy だけ離れた点を Y ,z 軸に沿っ
て Δz だけ離れた点を Z としたとき,点 P が点 P に移動すると点 X, Y , Z も,それぞれ点 X ,
Y , Z へと移動する.今,変形を二段階に分けて考える.第一段階では,点 X, Y, Z は座標軸に
平行である状態を保ったまま,点 X , Y , Z へと移動したと考える.第二段階では,点 X , Y , Z へと移動したと考える.第一段階は,点 P のまわりの点が,座標軸に沿って伸び縮みすることを
表し,垂直ひずみとして表される.第二段階は, xy 面,yz 面,zx 面に平行な面がゆがむことを
示しており,せん断ひずみとして表現される.以下,垂直ひずみ,せん断ひずみの順に定式化を
行う.
まず,図 2.6 の場合について,xy 面への投射を考える.図 2.7 に,xy 面へ投射した場合の模式
図を示す.
図 2.7: xy 面内の変形状態
9
変形前の線分 P X, P Y は,伸びや縮みによって P X , P Y へと移動する.今,線分 P X の
長さを Δx,P Y の長さを Δy ,P X の長さを Δx ,P Y の長さを Δy とする.P X の伸びは
Δx − Δx,P Y の伸びは Δy − Δy となる.変形の程度を表すひずみは,単位長さあたりの伸び
で表され,x 軸方向について式 (2.9),y 軸方向については式 (2.10) で示される.
εxx =
Δx − Δx
Δx
(2.9)
εyy =
Δy − Δy
Δy
(2.10)
また,図 2.7 に示すように,線分 P X と線分 P X が θx だけ傾いており,線分 P Y と線分
P Y が θy だけ傾いる場合,P X と P Y が作る xy 面に平行な面からのゆがみをせん断ひずみ
γxy として,θx + θy で表す.式 (2.11) に,せん断ひずみ γxy を示す.
γxy = θx + θy
(2.11)
一方,yz 面,zx 面についても同様に考えることで, z 軸方向の垂直ひずみ εzz ,yz 面のせん断
ひずみ γyz ,zx 面のせん断ひずみ γzx を得ることができる.
次に,ひずみと変位間の関係について述べる.簡単のため,x 軸方向についてのみ考える.点 P
の x 座標を Px とすると,点 X の x 座標は Px + Δx となる.また,点 P の x 座標を Px として,
式 (2.12) を満たす変位 ux (x, y, z) を定義する.
Px = Px + ux (x, y, z)
(2.12)
点 X の x 座標は,Px + x + ux (x + Δx, y, z) であるから,P X の長さ Δx は,変位ベクトル
の x 成分を用いて式 (2.13) で示される.
Δx = Δx + ux (x + Δx, y, z) − ux (x, y, z)
(2.13)
Δx を微小量として, 2 次以上の項を無視することで式 (2.14) を得る.
Δx ≈ Δx +
∂ux
Δx
∂x
(2.14)
式 (2.9) に式 (2.14) を代入することで,式 (2.15) を得ることができる.
εxx =
∂ux
Δx − Δx
=
Δx
∂x
(2.15)
同様に,垂直ひずみ εyy ,εzz についても,y 軸方向の変位 uy (x, y, z) と z 軸方向の変位 uz (x, y, z)
を用いて,それぞれ式 (2.16),式 (2.17) で定義できる.
εyy =
∂uy
Δy − Δy
=
Δy
∂y
(2.16)
εzz =
∂uz
Δz − Δz
=
Δz
∂z
(2.17)
次に,せん断ひずみ γxy について考える.今,点 P の y 座標は Py + uy (x, y, z) であり,点 X の y 座標は,Δx を微小量として二次以上の項を無視すると,式 (2.18) の関係が得られる.
Py + uy (x + Δx, y, z) ≈ Py + uy (x, y, z) +
10
∂uy
Δx
∂x
(2.18)
線分 P X と線分 P X の傾き θx は,式 (2.19) で示される.
θx =
Py + uy (x, y, z) +
∂uy
Δx − Py − uy (x, y, z)
∂uy
∂x
=
Δx
∂x
(2.19)
θy についても同様であり,式 (2.20) で表すことができる.
θy =
∂ux
∂y
(2.20)
式 (2.11) に,式 (2.19) と式 (2.20) を代入することで,せん断ひずみ γxy が得られる.式 (2.21)
に,せん断ひずみ γxy を示す.
γxy = θx + θy =
∂uy
∂ux
+
∂x
∂y
(2.21)
その他のせん断ひずみ γyz ,γz x についても,同様に得ることができる.式 (2.22) に γyz ,式
(2.23) に γzx を示す.
γyz =
∂uy
∂uz
+
∂y
∂z
(2.22)
γzx =
∂ux ∂uz
+
∂z
∂x
(2.23)
垂直ひずみとせん断ひずみを成分に持つベクトルであるひずみベクトル ε を導入し,ひずみと
変位の関係をベクトル表記すると式 (2.24) のようになる.
⎤
⎡
∂ux
⎥
⎢
∂x
⎥
⎢
⎥
⎢
∂uy
⎥
⎢
⎥
⎡
⎤ ⎢
⎥
⎢
∂y
εxx
⎥
⎢
⎥
⎢
⎥ ⎢
∂uz
⎥
⎢ εyy ⎥ ⎢
⎥
⎢
⎥ ⎢
∂z
⎥
⎢ εzz ⎥ ⎢
⎥
⎢
⎥
⎢
(2.24)
ε=⎢
= ⎢ ∂u
⎥
⎥
∂u
y
x
⎥
⎢ γxy ⎥ ⎢
+
⎢
⎥ ⎢
∂x ⎥
⎥
⎣ γyz ⎦ ⎢ ∂y
⎥
⎢ ∂u
∂u
⎥
⎢
y
z
γzx
+
⎥
⎢
⎢ ∂z
∂y ⎥
⎥
⎢
⎢ ∂uz
∂ux ⎥
⎦
⎣
+
∂x
∂z
式 (2.24) を,変位ーひずみ関係式と呼ぶ.
11
2.2.3
応力ーひずみ関係
線形弾性体では,応力とひずみには線形関係があり,式 (2.25) のように成分表示される.ただ
し,E はヤング率,ν はポアソン比である.ヤング率 E は,一方向の引張応力,または,圧縮応
力と応力方向に対するひずみの比であり,ポアソン比 ν は,荷重方向の伸びと荷重方向に直角方
向の縮みの比である.
⎤
⎡
ν
ν
⎡
⎤
1
0
0
0
1−ν 1−ν
⎥ εxx
⎢
⎥⎢
⎢ ν
⎥
ν
⎥⎢
⎢
⎥
1
0
0
0
⎥⎢
⎢
⎡
⎤
⎥ ⎢ εyy ⎥
⎢ 1−ν
1−ν
⎥
σxx
⎥⎢
⎢
⎥
ν
⎥⎢
⎢ ν
⎢
⎥
⎥
⎥⎢
⎢
⎢ σyy ⎥
1
0
0
0
⎥ ⎢ εzz ⎥
⎢ 1−ν 1−ν
⎢
⎥
⎥
⎥⎢
⎢
⎢ σzz ⎥
E(1 − ν)
⎥
⎥
⎢
⎢
⎥=
⎢
⎥
1 − 2ν
⎥
⎢ τ ⎥ (1 + ν)(1 − 2ν) ⎢ 0
⎢
0
0
0
0
⎥ ⎢ γxy ⎥
⎢
⎢ xy ⎥
⎥
2(1
−
ν)
⎥⎢
⎢
⎢
⎥
⎥
⎥⎢
⎢
⎣ τyz ⎦
⎥
⎥⎢
⎢
1 − 2ν
⎥⎢ γ ⎥
⎢ 0
0
0
0
0
τzx
⎥ ⎢ yz ⎥
⎢
2(1 − ν)
⎥
⎥⎣
⎢
⎦
⎥
⎢
1 − 2ν ⎦
⎣ 0
0
0
0
0
γzx
2(1 − ν)
(2.25)
2.3
破壊現象
本節では,破壊現象の分類と破壊発生の基準について記す.
2.3.1
脆性破壊と延性破壊
破壊現象は,物体内の原子間の結合が分離する現象であり,大きく分けて原子間の結合が引張
分離することに由来する脆性破壊 (Brittle fracture) と,せん断分離することに由来する延性破壊
(Ductile fracture) に大別される.図 2.8 に,破壊の分類を示す.
Object
Fracture
Tensile fracture
Cup and cone fracture
Shear fracture
Brittle fracture
Point fracture
Ductile fracture
図 2.8: 破壊の分類
12
脆性破壊は,最終破断までに著しい延びや絞りを伴うことのない破壊であり,巨視的にみると
引張破壊 (Tensile fracture) やカップコーン破壊 (Cup and cone fracture) が典型例である.また,
脆性破壊は分離破断とも言う.一方で,最終破断までに著しい伸びや絞りを伴う破壊を延性破壊
と呼ぶ.延性破壊の代表例としては,せん断破壊 (Shear fracture) や点状破壊 (Point fracture) が
挙げられる.
2.3.2
破壊基準
破壊は物体の応力が破断応力値 (Fracture stress) と呼ばれる臨界値を超えることにより生じる
が,破壊を引き起こす応力の違いによって三つの破壊条件,最大主応力条件 (Maximum principle
stress criterion),最大せん断応力条件 (Maximum shear stress criterion),せん断ひずみエネルギ
条件 (Shear strain energy criterion) がある.以下では,主応力を用いて各条件を記述するが,物
体のある点における最大主応力を σ1 ,最小主応力を σ3 ,中間主応力を σ2 とする.
• 最大主応力条件
破壊は主応力の蓄積によって引き起こされるという条件であり,一般に脆性破壊の破壊条
件として用いられる.破断応力値を σf p として,式 (2.26) の条件を満たす場合に破壊が生
じる.
σ1 ≥ σf p
(2.26)
最大主応力は,最大引張応力とも呼ばれ,最大主応力条件は,物体に対して引張操作を
行ったときの破壊基準として用いられる.また,破断面は,最大主応力の方向に垂直な面に
なると言われている.図 2.9 に,最大主応力条件による破壊の概念図を示す.
External force vector
Manipulated point
Object
Principal stress
max
Brittle fracture
min
Stress concentration
(Maximum principal stress exceed fracture stress)
図 2.9: 最大主応力条件による破壊の概念図
• 最大せん断応力条件
最大せん断応力条件は,経験的に延性破壊の破壊条件とされている.破断応力値を σf s と
して,式 (2.27) を満たす場合に破壊が生じるとされる.
σ1 − σ3 ≥ σf s
13
(2.27)
最大せん断応力条件は,物体に対しせん断操作を行ったときの破壊基準として用いられて
いる.一般に,せん断応力は,破断面に平行に作用すると言われている.図 2.10 に,最大せ
ん断応力条件による破壊の概念図を示す.
External force vector
Manipulated point
Object
Share stress
max
Ductile fracture
Stress concentration
min
(Maximum share stress exceed fracture stress)
図 2.10: 最大せん断応力条件による破壊の概念図
• せん断ひずみエネルギ条件
ミーゼス条件 (Mises criterion) とも言われ,せん断ひずみエネルギが破断応力値を超えた
場合に,破壊が発生するという破壊条件である.一般には延性破壊の破壊条件と言われてい
る.破断応力値を σf m として,式 (2.28) を満たす場合に,破壊が生じるとされる.
(σ1 − σ2 )2 + (σ2 − σ3 )2 + (σ3 − σ1 )2 ≥ 2σf2m
(2.28)
せん断ひずみエネルギ条件は,物体に対して引張操作とせん断操作を含むような操作を
行っているときに用いられる基準である.ただし,上述の二つの基準とは異なり,せん断ひ
ずみエネルギはスカラ値であり,破断面の方向は決定できない.図 2.11 に,せん断ひずみエ
ネルギ条件による破壊の概念図を示す.
External force vector
Manipulated point
Object
Mises stress
max
Ductile fracture
Stress concentration
min
(Maximum mises stress exceed fracture stress)
図 2.11: せん断ひずみエネルギ条件による破壊の概念図
14
2.3.3
破壊力学の知見
破壊現象は,破壊力学で定量的に扱われている.破壊力学では,破壊現象をき裂 (Crack) の発
生,および,進展現象であるとして,き裂先端近傍の応力やひずみを定式化している.簡単のた
めに,等方性で,き裂を有した厚みのある板状の線形弾性体について記す.図 2.12 に,説明に用
いるき裂を有する弾性板を示す.ただし,き裂の長さは 2a であり,き裂先端の曲率半径は限りな
く 0 に近いとする.また,厚み方向に関しては一様であるとする.
Object
Crack
2
図 2.12: 長さ 2a のき裂を有する弾性板
破壊力学では,代表的な三つの負荷条件下でき裂先端の応力場を定式化している.代表的な三
つの負荷条件というのは,いずれも,き裂の無限遠方で作用する一様な応力を仮定しており,き
裂面に垂直な方向に引張応力が作用している条件(モード I),き裂面に対して平行にせん断応力
が作用している条件(モード II),せん断応力が厚み方向に作用している条件(モード III)であ
る.図 2.13 に,三つの負荷条件を示す.
II
I
Crack
III
Crack
Crack
2
2
2
II
I
(a) (b) (c)
図 2.13: 代表的な三つの負荷条件
(a) モード I (b) モード II (c) モード III
15
III
モード I における引張応力を σI , モード II におけるせん断応力を τII , モード III におけるせん断
応力を τIII としたとき,各モードでの応力拡大係数 KI , KII , KIII というパラメータを定義する.
式 (2.29) に,応力拡大係数を示す.
√
KI = σI πa
√
KII = τII πa
√
KIII = τIII πa
(2.29)
破壊力学では,応力拡大係数を用いて,き裂先端近傍の応力を定式化している.ただし,先にも
記したように,厚み方向に関しては一様であるために,厚み方向に垂直な平面(x − y 平面)で考
える.図 2.14 に示すように,き裂先端を原点としてき裂の方向に軸をとった極座標系を考え,原
点からの距離を r ,偏角を θ とする.ただし,き裂先端は両端について定義されるが,簡単のため
に片側の先端についてのみ考える.
Object
Crack
2
図 2.14: き裂先端を原点とした極座標系
式 (2.30) にモード I の場合,式 (2.31) にモード II の場合,式 (2.32) にモード III の場合につい
て記す.
16
σx =
σy =
τxy =
σx =
σy =
τxy =
τyz =
τzx =
θ
θ
3θ
KI
√
cos
1 − sin sin
2
2
2
2πr
θ
θ
3θ
KI
√
cos
1 + sin sin
2
2
2
2πr
θ
KI
θ
3θ
√
cos sin cos
2
2
2
2πr
θ
KII
θ
3θ
sin
−√
2 + cos cos
2
2
2
2πr
θ
θ
3θ
K
√ II sin cos cos
2
2
2
2πr
θ
θ
3θ
K
√ II sin
1 − sin sin
2
2
2
2πr
θ
K
√ III cos
2
2πr
θ
KIII
sin
−√
2
2πr
(2.30)
(2.31)
(2.32)
式 (2.30)-(2.32) では,全ての応力が r の平方根の逆数に比例する形となっており,き裂先端に
近づくにつれて応力が大きくなることを示している.また,き裂先端においては,r = 0 となり応
力値が無限大に発散する点(特異点)として扱われる.
以上,破壊力学におけるき裂先端近傍の応力について示したが,無限遠の一様応力を仮定する
等,条件が固定されているために,任意の操作を行う場合に,破壊力学の知見をそのまま導入す
ることは不可能である.しかし,破壊現象を,き裂の発生,および,進展現象として扱うことや
き裂の先端では応力が大きくなるといったことは,破壊現象を取り扱う上では重要な事項である.
17
第 3 章 力覚提示技術と応用
本章では,力覚提示可能な VR システムの要素技術と応用例について記す.まず,物体操作を
対象とした多くの力覚提示 VR システムで用いられている有限要素モデル [38]-[42] について記す.
次に,ユーザへの反力提示を可能にするために必要な,力覚提示デバイスについて記す.最後に,
既存研究において,破壊現象を含む操作を可能とする研究,および,記録・再生型の反力提示を
行う研究について例を示して紹介する.
3.1
有限要素モデル
本節では,等方性の線形弾性体を対象とする有限要素法の説明,および,実時間で有限要素モ
デルの操作を可能とする手法について記す.
3.1.1
支配方程式
外力の作用する弾性体について,以下に示す仮想仕事の原理が成立する.
「外力作用の元で釣り合い状態にある弾性体の各点に,任意の微小な仮想変位を与えたとき,外
力と内力が仮想変位によって行われる仕事の総和は 0 である.
」
静的な力の釣り合いの元では,仮想仕事の原理により,各時刻に外力と内力の釣り合いを仮定
し,各時刻に成立する剛性方程式を解くことで,変位と反力を求める.
外力は表面力と体積力に分けられる.表面力とは,物体が他の物体と接触して物体表面に及ぼ
す力であり,例えば,水中にある物体の表面に作用する水圧などが表面力に該当する.また,体
積力とは,物体そのものに対して,物体の体積に比例して作用する力であり,例としては重力が
挙げられる.
以下では,体積力から受ける影響を考慮せずに表面力のみを考える.単位面積当りの表面力ベ
クトル F を式 (3.1),仮想変位ベクトル u∗ を式 (3.2),仮想ひずみベクトル ε∗ を式 (3.3) で定義
する.
⎡
⎤
Fx
⎢
⎥
F = ⎣ Fy ⎦
Fz
(3.1)
⎡
⎤
u∗x
⎢
⎥
u∗ = ⎣ u∗y ⎦
u∗z
18
(3.2)
⎡
⎢
⎢
⎢
⎢
∗
ε =⎢
⎢
⎢
⎢
⎣
ε∗xx
ε∗yy
ε∗zz
∗
γxy
∗
γyz
∗
γzx
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(3.3)
式 (3.1)-(3.3) と,応力ベクトル σ を用いて,仮想仕事の原理を定式化すると式 (3.4) のように
なる.なお,右上付き文字の T は転置を表す.
(3.4)
u∗ T F dS − ε∗ T σdV = 0
式 (3.4) で示される仮想仕事の原理式と式 (2.24) の変位ーひずみ関係式, および式 (2.25) の応
力ーひずみ関係式の三式が有限要素法の支配方程式となる.
3.1.2
剛性方程式導出
三つの支配方程式から,剛性方程式を導出する過程を四面体要素を用いる場合について説明す
る.図 3.1 に,有限要素法の概念図を示す.なお,図 3.1 においては,簡単のため,二次元的に物
体を表現し,四面体要素は三角形要素としている.
Object
Finite element model
Node
Element
図 3.1: 有限要素法の概念図
19
離散化するにあたり以下の二つの仮定を設ける.
1. 応力,ひずみは要素内で一定である.
2. 外力や変位は節点を介してのみ伝わる.
まず,1 番の仮定は要素内の変位が線形に変化することと等価である.つまり,四面体内の変位
場を式 (3.5) のように仮定することと等しい.
ux = α0 + α1 x + α2 y + α3 z
uy = α4 + α5 x + α6 y + α7 z
(3.5)
uz = α8 + α9 x + α10 y + α11 z
ただし,ux は x 方向の変位, uy は y 方向の変位, uz は z 方向の変位である.また,α0 , α1 , · · · , α11
は定数である.四つの頂点を i, j, k, l,頂点変位を (uxi , uyi , uzi ), (uxj , uyj , uzj ), (uxk , uyk , uzk ),
(uxl , uyl , uzl ),頂点座標を Vi (xi , yi , zi ), Vj (xj , yj , zj ), Vk (xk , yk , zk ), Vl (xl , yl , zl ) として式 (3.5) に
代入すると,x 方向の変位について,式 (3.6) が得られる.
⎡
⎤ ⎡
⎤⎡
⎤
uxi
1 xi yi zi
α0
⎢ u ⎥ ⎢ 1 x y z ⎥⎢ α ⎥
⎢ xj ⎥ ⎢
j
j
j ⎥⎢ 1 ⎥
(3.6)
⎢
⎥=⎢
⎥⎢
⎥
⎣ uxk ⎦ ⎣ 1 xk yk zk ⎦ ⎣ α2 ⎦
uxl
1 xl yl zl
α3
未知数係数 α0 , α1 , α2 , α3 についてみると,式 (3.7) が成立する.
⎡
⎡
⎤ ⎡
⎤−1 ⎡
⎤
ai aj ak
α0
1 xi yi zi
uxi
⎢
⎢ α ⎥ ⎢ 1 x y z ⎥ ⎢ u ⎥
1 ⎢ bi bj bk
⎢ 1 ⎥ ⎢
⎢ xj ⎥
j
j
j ⎥
⎢
⎢
⎥=⎢
⎥ ⎢
⎥=
⎣ α2 ⎦ ⎣ 1 xk yk zk ⎦ ⎣ uxk ⎦ 6V ⎣ ci cj ck
α3
1 xl yl zl
uxl
di dj dk
al
bl
cl
dl
⎤⎡
⎥⎢
⎥⎢
⎥⎢
⎦⎣
uxi
uxj
uxk
uxl
⎤
⎥
⎥
⎥
⎦
(3.7)
ただし,V は四面体の体積である.式 (3.7) から求まる α0 , α1 , α2 , α3 を式 (3.5) に代入すると,
式 (3.8) が得られる.
ux =
1
{(ai + bi x + ci y + di z)uxi + (aj + bj x + cj y + dj z)uxj
6V
+(ak + bk x + ck y + dk z)uxk + (al + bl x + cl y + dl z)uxl }
(3.8)
同様に y 方向について式 (3.9),z 方向について式 (3.10) が成り立つ.
uy =
uz =
1
{(ai + bi x + ci y + di z)uyi + (aj + bj x + cj y + dj z)uyj
6V
+(ak + bk x + ck y + dk z)uyk + (al + bl x + cl y + dl z)uyl }
(3.9)
1
{(ai + bi x + ci y + di z)uzi + (aj + bj x + cj y + dj z)uzj
6V
+(ak + bk x + ck y + dk z)uzk + (al + bl x + cl y + dl z)uzl }
(3.10)
また,式 (3.8), 式 (3.9),および,式 (3.10) において,節点変位成分の係数部分を式 (3.11) のよ
うに置き換える.
Ni = ai + bi x + ci y + di z
Nj = aj + bj x + cj y + dj z
Nk = ak + bk x + ck y + dk z
Nl = al + bl x + cl y + dl z
20
(3.11)
結果的に,式 (3.12) が得られる.
ux = Ni uxi + Nj uxj + Nk uxk + Nl uxl
uy = Ni uyi + Nj uyj + Nk uyk + Nl uyl
(3.12)
uz = Ni uzi + Nj uzj + Nk uzk + Nl uzl
変位ーひずみ関係式 (2.24) に式 (3.12) を代入すると,式 (3.13) が得られる.
⎡
⎤
⎡
∂Ni
∂Nj
∂Nk
∂Nl
0
0
0
0
0
0
0 ⎥⎢
⎢
⎢ ∂x 0
∂x
∂x
∂x
⎥⎢
⎢
⎥
⎢
⎢
∂Nj
∂Nk
∂Nl
⎥⎢
⎢ 0 ∂Ni 0
0
0
0
0
0
0
⎥⎢
⎢
∂y
∂y
∂y
∂y
⎥⎢
⎢
⎥⎢
⎢
⎢
∂Nj
∂Nk
∂Nl ⎥ ⎢
∂Ni
⎥⎢
⎢ 0
0
0
0
0
0
0
0
⎢
⎢
∂z
∂z
∂z
∂z ⎥
⎥⎢
⎢
[ε] = ⎢
⎥
⎢
⎢ ∂Ni ∂Ni 0 ∂Nj ∂Nj 0 ∂Nk ∂Nk 0 ∂Nl ∂Nl 0 ⎥ ⎢
⎥
⎢
⎢ ∂x ∂y
∂x ∂y
∂x ∂y
∂x ∂y
⎥⎢
⎢
⎥⎢
⎢
⎢
∂Nj ∂Nj
∂Nk ∂Nk
∂Nl ∂Nl ⎥ ⎢
∂Ni ∂Ni
⎥⎢
⎢ 0
0
0
0
⎢
⎢
∂y ∂z
∂y ∂z
∂y ∂z
∂y ∂z ⎥
⎥⎢
⎢
⎦
⎣ ∂Ni
∂Nj ∂Nk
∂Ni ∂Nj
∂Nk ∂Nl
∂Nl ⎢
⎣
0
0
0
0
∂x
∂z ∂x
∂z ∂x
∂z ∂x
∂z
uxi
uyi
uzi
uxj
uyj
uzj
uxk
uyk
uzk
uxl
uyl
uzl
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
(3.13)
式 (3.13) の節点変位ベクトルの係数行列を B として,四面体の節点変位ーひずみ関係式は,節
点変位ベクトル un を用いて,式 (3.14) のように書くことができる.なお,以下では,節点に関す
る量には添字 n を用いる.
ε = Bun
(3.14)
次に,2 番目の仮定から要素に作用する表面力 F と同等の仕事をする節点力を定義する.つま
り,節点力ベクトル f n と節点変位ベクトル un について,式 (3.15) が成立するように設定する.
(3.15)
u∗ T F dS = u∗n T f n
また,要素内で応力,ひずみが一定であることから式 (3.16) が成立する.
ε∗ T σdV = V ε∗ T σ
(3.16)
式 (3.15) と式 (3.16) から仮想仕事の原理式は式 (3.17) のようになる.
u∗n T f n = V ε∗T σ
(3.17)
また,変位ーひずみ関係式と応力ーひずみ関係式について,仮想ひずみベクトル ε∗ と,仮想変
位ベクトル u∗n の関係は式 (3.18),応力ベクトル σ と,節点変位ベクトル un の関係は式 (3.19) で
示される.
ε∗ = Bu∗n
(3.18)
σ = Dε = DBun
(3.19)
21
式 (3.17) に式 (3.18) と式 (3.19) を代入すると式 (3.20) が得られる.
u∗n T f n = V ε∗T σ
= V (Bu∗n )T DBun
= u∗n T (V B T DB)un
(3.20)
任意の仮想節点変位ベクトルの転置ベクトル u∗n T について,式 (3.20) が成立するためには式
(3.21) を満たす必要がある.
f n = (V B T DB)un
(3.21)
K e = V B T DB とおくと,最終的に要素剛性方程式 (3.22) が得られる.
f n = K e un
(3.22)
また,物体を構成する全要素剛性方程式をまとめると式 (3.23) の剛性方程式が導出される.
f = Ku
3.1.3
(3.23)
実時間インタラクション
次に,有限要素法の計算量について述べる.まず,式 (3.23) の物体の剛性方程式を,固定節点
(Fixed node : 添字 i),自由節点 (Free node : 添字 r),接触節点 (Contact node : 添字 c) につい
ての部分ベクトル,もしくは,部分行列を用いて示す.また,図 3.2 に,固定節点,自由節点,お
よび,接触節点の概念図を示す.
Contact node
Fixed nodes
Free nodes
図 3.2: 固定節点,自由節点,接触節点の概念図
22
成分表示すると式 (3.24) が得られる.
⎡
⎤ ⎡
fi
K ii K ir
⎢
⎥ ⎢
=
f
K
⎣ r ⎦ ⎣ ri K rr
fc
K ci K cr
⎤⎡
⎤
ui
K ic
⎥⎢
⎥
K rc ⎦ ⎣ ur ⎦
K cc
uc
(3.24)
初めに,静的な力の釣り合いから固定節点変位と外力,および,接触点以外の外力は 0 である
とし,式 (3.25) を得る.
0
K rr K rc
ur
=
(3.25)
K cr K cc
uc
fc
式 (3.25) において,既知数(入力)は接触点変位 uc であり,未知数(出力)は接触点外力 f c
と接触点以外の自由節点変位 ur である.
式 (3.26) と式 (3.27) に,式 (3.25) を連立方程式の形で示す.
0 = K rr ur + K rc uc
(3.26)
f c = K cr ur + K cc uc
(3.27)
全体の節点数を N , 接触節点数を Nc , 接触点以外の自由節点数を Nr , 固定節点数を Ni としたと
き,式 (3.26) は,3×(N −Ni −Nc ) 次の連立方程式となり,式 (3.27) は,次数 (3×Nc , 3×(N −Ni ))
の行列と次数 3 × Nc のベクトルの乗算となる.仮に,接触点数を一点であるとした場合,毎回の
インタラクション時に解く必要がある連立方程式の次数は 3 × (N − Ni − 1) となり,単純なモデ
ルでない限りは実時間計算が困難である.
しかし,式 (3.25) の剛性方程式を剛性マトリクスの逆行列を前処理段階で求めておくことで,毎
回の有限要素法の計算量を削減する手法が広田らによって提案されている [11].広田らの手法で
は,式 (3.28) に示す剛性マトリクスの逆行列 L を用いて,式 (3.25) を式 (3.29) のように変形する.
L = K −1
ur
uc
=
Lrr Lrc
Lcr Lcc
(3.28)
0
fc
(3.29)
連立方程式の形にすると,式 (3.30) と式 (3.31) が得られる.
uc = Lcc f c
(3.30)
ur = Lrc f c
(3.31)
式 (3.30) は,3 × Nc 次の連立方程式であり,式 (3.31) は,次数 (3 × Nc , 3 × Nr ) の行列と次数
3 × Nc のベクトルの乗算となる.先と同じように,接触点数を一点と仮定した場合には,毎回の
有限要素法の計算で三次の連立方程式を解くだけでよいことになり,逆行列を用いない場合に比
べて計算量が少ない.
上記では,節点に接触することを前提として,話を進めてきたが,実際の操作では,面に接触
することもあり得る.面に接触した場合は,面を構成する節点に独立に変位を与え,式 (3.29) で
23
Polygon P
Reaction force vector
of point
Displacement vector
Reaction force vector
of vertex
Displacement vector
Reaction force vector
of vertex
Reaction force vector
of vertex
図 3.3: 面上の任意の点に接触した場合の反力計算
反力を計算したのち,得られた反力を面積座標で重み付けして補間することで,面上の任意の点
における反力を算出する.図 3.3 に,面上の任意の点に接触した時における,反力生成の模式図を
示す.
今,節点 a, b, c で構成される面 P 上の点 q に接触し,ベクトル uq だけ変位を与えた場合を考
える.節点 a, b, c に対して,独立にベクトル uq だけ変位を与えた場合の反力ベクトルを f a , f b ,
f c とするとき,接触点 q における反力 f q は,上記三つの反力ベクトルに,対応する面積座標で重
み付け和で与えられる.面積座標とは,三角形内の一点が作る三つの三角形の面積の比率である.
図 3.4 に,節点 a, b, c に対応する面積座標を Aa , Ab , Ac とした場合の例を示す.
図 3.4: 面積座標の例
面積座標間には式 (3.32) が成立する.
Aa + Ab + Ac = 1
24
(3.32)
接触点 q における反力ベクトル f q は,式 (3.33) で与えられる.
f q = Aa f a + Ab f b + Ac f c
3.2
(3.33)
力覚提示デバイス
本節では,代表的な力覚提示デバイスについて特徴をまとめる.
3.2.1
接地型
接地型の力覚提示デバイスの例として,PHANToM, Falcon, SPIDAR, HapticMaster の四つを
示す.本小節で紹介する接地型力覚提示デバイスは,ペン型や球形の把持部を介してユーザに力
覚を提示するデバイスである.一方で.デスクトップ型のデバイスであるために,使用場所があ
る程度限られることや一点接触を想定して作られており,複数点接触を想定するアプリケーショ
ンなどでは複数台のデバイスが必要になるといった欠点がある.
PHANToM
SenseAble 社製 PHANToM[43] は,アーム (Arm) の多関節機構を駆動機構とする多自由度の力
覚を提示可能なデバイスである.ペン (Pen) 型の把持部を介してユーザに力覚を提示する.図 3.5
に三自由度の力覚提示が可能な PHANToM Omni の概観,図 3.6 に六自由度の力覚提示可能な
PHANToM Premium の概観を示す.
Pen
Arm
Arm
図 3.5: PHANToM Omni の概観 [43]
Pen
図 3.6: PHANToM Premium の概観 [43]
PHANToM はペン形状の把持部を,仮想的な道具や人間の指とするようなアプリケーションに
用いられ,現在最も普及している力覚提示デバイスである.
Falcon
Novint 社製 Falcon[44] は,胴体部と把持部 (Grip) が三本のアームで連結されている力覚提示デ
バイスである.入出力ともに三自由度であり,球形の把持部を介して力覚を提示する.図 3.7 に,
Falcon の外観を示す.
25
Arm
Grip
図 3.7: Falcon の概観 [44]
Falcon はゲーム用のコントローラとして開発されており,前述の PHANToM に比べて安価では
あるが,可動域が小さく,広い可動域が必要される用途には不向きである.
SPIDAR
SPIDAR[45] は佐藤らによって開発された力覚提示デバイスである.図 3.8 に外観を示す.
String
Grip
図 3.8: SPIDAR の概観 [45]
SPIDAR は立方体のフレームの各隅から中央の把持部へ糸 (String) が張られており,モータで
糸の張力を制御することで,位置入力と反力提示を可能とするデバイスである.六自由度の力覚
提示が可能である.
SPIDAR は,力覚提示部が糸で構成されているため,遮蔽が少なく,視覚空間と力覚空間の統
合が行いやすいことや,両手を用いて操作することができるため,現実感のある操作が可能であ
るといった利点がある.一方で,操作によっては手や指とワイヤが干渉する可能性があることや,
作業領域内に手や指以外の物,例えばディスプレイや頭部などを入れることができないといった
欠点もある.
26
HapticMaster
HapticMater[46] は六自由度の入力が可能な力覚提示デバイスである.図 3.9 に外観を示す.
Pantograph
Grip
図 3.9: HapticMaster の概観 [46]
HapticMaster は,三つのパンタグラフ (Pantograph) それぞれに三つの DC モータが取り付け
られ,パルス幅変調方式によりトルクの制御を行っている.比較的広い可動域を持つが,デバイ
ス自体が大型であるため,使用場所が限られる欠点がある.
3.2.2
装着型
CyberGrasp・CyberForce
Immersion 社で開発された CyberGrasp[47] は,グローブ状の装着型力覚提示デバイスである.
また,CyberForce は,ロボットアーム (Robot arm) の先端に CyberGrasp を取り付けたデバイス
である.CyberGrasp の概観を図 3.10(a),CyberForce の概観を図 3.10(b) に示す.
Robot arm
Wire
Cyber grasp
(a) (b)
図 3.10: CyberGrasp と CyberForce の概観 [47]
(a)CyberGrasp (b)CyberForce
27
CyberGrasp は,各指に力覚提示装置を装着し,ワイヤ (Wire) 駆動により独立した力覚提示を
行うデバイスである.CyberGrasp では,仮想物体に触れる,把持する感覚が得られることや,作
業空間が広いことが利点である.また,CyberForce は,手の動きを三次元の位置と角度の変化と
してとらえ,手の甲に力を提示することにより,例えば,仮想空間内で自動車のハンドルを握った
ときの重量感といった感覚を提示できる.接地型の力覚提示デバイスに比べて,握る動作に関し
てより現実感のある力覚提示が可能ではあるが,装着感を取り除くことはできない.
Gravity Grabber
Gravity Grabber[48] は,仮想の箱やグラスの中に物体が入った時の重さを体感できるデバイス
である.Gravity Grabber の概観を,図 3.11 に示す.
Motor
図 3.11: Gravity Grabber[48]
Gravity Grabber では,二個のモータ (Motor) でベルトを駆動することにより反力を生成する.
二個のモータが同方向に回転すると,指とベルトの接触面に水平な方向の力を提示可能であり,二
個のモータが逆方向に回転すると,指とベルトの接触面に垂直な方向の力を提示可能である.
流体グローブ
流体グローブ [49] は,空気圧ベローズ (Bellows) を用いて力を生成する装着型の力覚提示デバイ
スである.図 3.12 に,流体グローブの概観を示す.
Bellows
図 3.12: 流体グローブ [49]
28
流体グローブは,空気圧ベローズ内部に空気を送り込み,加圧して剛性を増すことで,指に拘
束力を与えるという動作原理である.流体グローブでは,小指以外の四指に空気圧ベローズを取
り付けることで,各指に独立した力覚を提示し,物体を把持する力覚を提示できる.
本節では,種々の力覚提示デバイスを紹介した.本研究では,固体の仮想物体に対して,ユー
ザは一点を押す,もしくは一点を掴み引っ張るといった操作を仮定する.一点の入出力と,三自由
度の入出力が要求される.要件を満たすデバイスの中で,球体の把持部を掴む操作,および,把
持部に存在するボタンを押すといった操作が,実際に物体を掴む操作を連想させることができる
と考え,Novint 社 Falcon を用いる.
3.3
アプリケーション
本節では,対話的に仮想物体の破壊現象を扱っているモデル,および,記録・再生型の反力提
示を行っているモデルについて例を示して説明する.
3.3.1
対話的破壊を扱うモデル
仮想切断シミュレータ
仮想空間内で,ナイフなどの切断器具を用いて物体を切断する操作を模擬できるシミュレータ
である [50]-[54].図 3.13 に,仮想空間内の膜状の物体 (Virtual membrane) を仮想ナイフ (Virtual
knife) で切断している様子を示す.
Virtual knife
Virtual membrane
図 3.13: 仮想切断シミュレータ [50]
仮想物体を,質点とバネ,および,ダンパで結合されたモデルとして構築しており,粘弾性体
の破壊現象を実時間で模擬可能である.しかし,応力やひずみなどの物理量の計算や,ヤング率
やポアソン比などの材料定数の適用が困難であるという欠点がある.また,仮想切断器具と物体
との接触点で破壊が発生し,切断器具の進む方向に破壊が進行するという制約を設けているため
に,厳密な物理モデルとは言えない.
剥離手技シミュレータ
仮想物体の応力状態に基づいた破壊を可能にする研究例として,粂らが開発した生体軟組織剥
離手技シミュレータ [55] がある.図 3.14 に,剥離シミュレータの概観を示す.
29
図 3.14: 剥離シミュレータ [55]
剥離手技のシミュレータでは,物体の応力状態に基づいた破壊を実現している.しかし,幾何
形状や剛性マトリクスの再構築計算に高速化の手法を用いることなく,物体全体の再構築計算を
行っているために,単純なモデルでない限り実時間計算は困難である.
3.3.2
記録・再生型モデル
Haptic texture
Haptic texture[56] は,物体の場所ごとに弾性値を割り当てる概念である.図 3.15 に,人体の表
面モデル (Human body surface) に弾性値 (Value of elasticity) を割り当てた例を示す.弾性値が
高いほど白く,低いほど灰色になっている.
Human body surface
Value of elasticity
high
low
図 3.15: Haptic texture の概念図 [56]
30
Haptic texture では,前処理段階で,物体の表面の各点に変位を与えたときの反力値から,各
点における弾性値を算出し割り当てる.実時間処理中には,操作に応じて割り当てられた弾性値
を元に反力を算出する.物理計算を行うことなく,提示反力が計算できるために,実時間処理中
に行う物理計算のコストを削減することが可能である.
一方で,Haptic texture の概念だけでは,操作による物体の変形について,何も規定していな
いために,対話的な物体操作を仮定するシステム構築には不十分である.しかし,前処理段階で
パラメータを記録し,実時間処理時に再生することで,実時間計算をする物理計算量を削減する
という概念は,以下に紹介するアプリケーションや本研究においても重要な概念となっている.
インパルス応答の記録・再生モデル
田川らの研究 [57], [58] では,動的な変形モデルに対して,記録・再生型のアプローチを適用し
ている.前処理段階で,表面節点にインパルス的な変位を与えたときの変形と反力の応答を記録
し,実時間処理時には,入力変位と記録した応答を畳み込むことで,変形と反力を計算している.
図 3.16 に,田川らのモデルに変位 (Displacement) を与えたときの,変形の様子を示す.
Displacement
Time
図 3.16: 記録・再生型アプローチを用いた変形モデル [58]
田川らの研究では,変形操作に対象を限定しており,破壊現象を取り扱うことはできない.し
かし,前処理段階における処理を工夫することで,シミュレーション中に行う物理計算量を削減
することができるという重要な事実を示す研究である.
Haptic video
Haptic video[59] とは,熟練者の作業を記録し,訓練者に力覚提示デバイスを用いて提示するこ
とで技術指導を行うシステムである.従来,技術指導は知識の伝達が主な手段であったが,Haptic
Video によって,身体運動を伴った効果的な技術指導が可能となる.また,熟練者の作業を記録す
るために,訓練者は熟練者のいない環境でも技術訓練が可能となる.図 3.17 に,Haptic video の
概念図を示す.
31
図 3.17: Haptic video の概念図 [59]
Haptic video では,熟練者 (Expert) が作業している時の動きや反力を,力覚提示デバイス (Haptic
device) を用いて記録するだけではなく,作業の様子をカメラ (camera) で撮影し,再生時には,訓
練者 (Trainee) の主体的操作による反力提示が行われるとともに, LCD により撮影された作業風
景も再生している.訓練者は,力覚のみならず視覚的にも,熟練者の作業を学ぶことが可能であ
る.
本研究でも,変形や破壊進行計算については実時間処理中に計算を行うが,実時間計算が困難
な物体全体のトポロジ更新や反力提示については記録・再生型のアプローチを採用する.
32
第 4 章 実時間トポロジ変化モデル
本章では,本研究で提案する対話的に破壊操作可能な有限要素モデルの構築手法について述べ
る.提案手法は実行順序で大別すると,前処理段階での剛性マトリクスの逆行列,および,反力
データの生成と,実時間処理時の物理計算,および,前処理データの補正の二つに分けられる.本
章では, 提案手法の全体像を記した後,両者について詳細を記す.
4.1
提案手法概要
本節では,提案手法の概要と処理の流れについて述べる.なお,仮想物体として,初期き裂を
有する等方性の線形弾性体を仮定する.まず,提案手法で用いる物理モデルについて記す.図 4.1
に提案手法で用いるモデルの概略図を示す.以後,簡単のため,提案手法の説明は,図 4.1 に示す
二次元のモデルを用いて行う.
Local object
Boundary nodes
Global object
Crack
図 4.1: 提案手法で用いるモデルの概略図
破壊進行時の再構築計算の計算コストを削減するために,再構築計算する要素数を削減する.提
案モデルでは,物体は物体全体から構成されるグローバルオブジェクト (Global object) と,き裂
(Crack) 近傍の局所領域から構成されるローカルオブジェクト (Local object) で構成されるとする.
ローカルオブジェクトとグローバルオブジェクトは,境界節点 (Boundary node) を共有しており,
毎回の再構築計算はローカルオブジェクトについてのみ行う.
33
Pre-calculation phase
Calculation of inverse stiffness matrices
Calculation of reaction forces
Real-time calculation phase
Images of an object
Rendering
User
Reaction force
Displacement
Reaction forces database
Global object
Displacement of
boundary nodes
Local object
s ma
ffnes
se sti
Inver
FEM1
FEM2
trix
Stress
Pre-calculated reaction forces data
Inverse stiffness matrices database
Determination of
crack propagation
Mesh and matrix restructuring
(Calculation of stiffness matrix)
Pre-calculated
Inverse stiffness matrices
FEM1 : Solves a stiffness equation with inverse stiffness matrix
FEM2 : Solves a stiffness equation as simultaneous equation
図 4.2: 処理の流れ
以下,処理の流れについて記す.図 4.2 に処理の流れを示す.
前処理段階 (Pre-calculation phase) では,多様なトポロジにおける物体全体の剛性マトリクス
の逆行列,および,多様なトポロジと複数の操作パターンにおける反力値を記録する.以下に,前
処理計算を行うために必要なパラメータを示す.また,括弧内の文字は,各パラメータのパター
ン数を表す.
34
• 破壊進展回数 (Ncp )
• 破壊進行角度 (Ncd )
• 破壊進行量 (Ncl )
• ユーザの操作点位置 (Nmp )
• ユーザの操作方向 (Nmd )
• ユーザの操作量 (Nml )
物体のトポロジは,破壊が進行した回数と破壊が進行した角度,および,破壊の進行量によっ
て決まるために,上記パラメータの中で,破壊進展回数と破壊進行角度,破壊進行量は剛性マト
リクスの逆行列と反力データについて共通のパラメータである.また,ユーザの操作は,ユーザ
の操作点の位置と操作した方向,さらに,操作した量によって定義できるために,反力限定のパ
ラメータとなる.生成されるデータ量は,物体全体の剛性マトリクスの容量を M (byte),反力の
自由度を 3 とし,倍精度 (8byte) を仮定すると,剛性マトリクスの逆行列の容量は式 (4.1),反力
データの容量は式 (4.2) で示される.
M×
Ncp
(Ncd × Ncl )i
(4.1)
i=1
3 × 8 × Nml × Nmp × Nmd ×
Ncp
(Ncd × Ncl )i
(4.2)
i=1
実時間処理中 (Real-time calculation phase) の処理については,まず,ユーザ (User) は力覚提示
デバイスを介して,グローバルオブジェクトに対して変位を与える.変位を入力として,グローバ
ルオブジェクトについて有限要素法計算 (FEM1) を行い,グローバルオブジェクトの全ての節点の
変位を計算する.次に,グローバルオブジェクトとローカルオブジェクトとの境界節点の変位を入
力として,ローカルオブジェクトに対して有限要素法計算 (FEM2) を行い,ローカルオブジェクト
の全ての節点変位を求める.得られた変位を元に,ローカルオブジェクト内で応力計算を行い,破
壊判定 (Determination of crack propagation) を行う.破壊判定の結果,破壊が進行するとなれば,
ローカルオブジェクトの剛性マトリクスや幾何形状の再構築計算 (Mesh and matrix restructuring)
を行い,ローカルオブジェクトを新しいき裂先端の領域近傍で再確保する.再確保の際,グロー
バルオブジェクトにローカルオブジェクトの再構築計算の結果を反映させる必要がある.提案手
法では,ローカルオブジェクトの再構築計算の結果に基づいて,前処理で記録した剛性マトリク
スの逆行列データ (Pre-calculated inverse stiffness matrices) の中から,一番類似するデータをグ
ローバルオブジェクトに適用する.しかし,類似の逆行列を適用するだけであり,実時間処理中
のトポロジを正確に考慮しているとは言えず,物体トポロジを詳細に反映した反力提示は困難で
ある.実時間処理中の物体トポロジとユーザ操作を元に,前処理した反力データ (Pre-calculated
reaction forces data) を補正してユーザに提示することで,物体トポロジを詳細に反映した反力提
示を行う.
35
4.2
再構築計算時間削減
提案モデルでは,破壊進行に合わせて,ローカルオブジェクトを更新する必要がある.破壊進
行は,ローカルオブジェクトの応力計算の結果に依存するため,ローカルオブジェクトの更新は
応力計算の結果から決定されることになる.
本節では,ローカルオブジェクトを用いた有限要素法計算と破壊進行を決定する方法,ならび
にローカルオブジェクトの更新方法について記す.また,破壊によるトポロジ変化に伴い生じる
幾何形状と剛性マトリクスの再構築計算についても述べる.
4.2.1
ローカルオブジェクトを用いた有限要素法計算
まず,ユーザからの変位入力を受けて,グローバルオブジェクトについて有限要素法計算を行
う.グローバルオブジェクトの剛性マトリクスの逆行列を L,操作点の添字を c, 境界節点の添字
を b,境界節点以外の自由節点の添字を o,u を変位ベクトル,f を外力ベクトルとすると,グロー
バルオブジェクトの剛性方程式は,式 (4.3) で示される.
⎡
⎤ ⎡
⎤⎡
⎤
uc
Lcc Lbc Loc
fc
⎢
⎥ ⎢
⎥⎢
⎥
⎣ ub ⎦ = ⎣ Lcb Lbb Lob ⎦ ⎣ f b ⎦
uo
Lco Lbo Loo
fo
(4.3)
式 (4.3) において,入力となる操作点変位ベクトル uc と,静的な力の釣り合いのため導入される
操作点以外の外力ベクトル f b , f o は零ベクトルであり既知である.計算順序は,はじめに式 (4.4)
により,操作点の外力ベクトル f c を計算する.
f c = Lcc −1 uc
(4.4)
次に,式 (4.5) により,操作点以外の自由節点の変位を求める.
ub
Lcb
=
fc
uo
Lco
(4.5)
続いてローカルオブジェクトでは,グローバルオブジェクトにおける有限要素法計算によって算
出された境界節点変位ベクトル ub を入力として,有限要素法計算が行われる.ローカルオブジェ
クトについて,成立する剛性方程式を式 (4.6) に示す.なお,境界節点を添字 b,その他の自由節
点を添字 o で表し,u は変位,f は外力,K は剛性マトリクス成分を示している.
fb
fo
=
K bb K ob
K bo K oo
ub
uo
(4.6)
式 (4.6) において,入力となる境界節点変位ベクトル ub と,静的な力の釣り合いを仮定するた
めに導入される自由節点外力 f o は零ベクトルは既知である.また,境界節点外力ベクトル f b と
自由節点変位ベクトル uo は未知である.ローカルオブジェクトにおける有限要素法計算は,ロー
カルオブジェクトの変位,および,応力を計算することが目的であるために,未知数の中で境界
節点外力ベクトル f b は計算する必要がない.結果として,解くべき方程式を式 (4.7) に示す.
0 = K bo ub + K oo uo
36
(4.7)
境界節点数を Nb ,自由節点数を No としたとき,式 (4.7) は 3 × No 次の連立方程式となる.
以下,計算時間について考える.前章では,有限要素法の高速計算手法として,剛性マトリク
スの逆行列を前処理段階で求める手法を紹介した.ローカルオブジェクトに対する有限要素法計
算においても,逆行列を前処理段階で求めることで,毎回の有限要素法計算のコストを削減でき
ると考えられる.
しかし,接触点変位を入力とする,グローバルオブジェクトにおける有限要素法の場合と異な
り,ローカルオブジェクトにおける有限要素法計算では,入力節点数は境界節点数となり,ロー
カルオブジェクトの全体節点数の中でも大きな割合を占めることとなる.従って,逆行列を前処
理段階で求める手法による高速化がグローバルオブジェクトの場合ほど期待できない.また,逆
行列を前処理で求める手法をローカルオブジェクトに採用した場合,再構築計算の度,ローカル
オブジェクト全体の剛性マトリクスの逆行列を計算する必要が生じる.一方,本節で紹介した方
法では,剛性マトリクスの逆行列計算を必要とはせず,剛性マトリクスの部分マトリクス K oo の
LU 分解のみを前処理することで実行できる. LU 分解計算は,実時間処理時に連立方程式を解く
時の計算量を削減するために行う.
以上のことから,本研究では,ローカルオブジェクトにおける有限要素法計算について,逆行
列計算を行わない.再構築計算時に行うべき処理は,剛性マトリクスの部分行列 K oo の LU 分解
計算となる.
4.2.2
ローカルオブジェクトの破壊進行計算
操作によって破壊が進行する場合を考える.まず,ユーザからの変位入力を受けて,ローカル
オブジェクト内で応力計算が行われた状態から説明を始める.図 4.3 に,ローカルオブジェクト内
で,有限要素法計算が行われた様子を模式的に二次元モデルで示す.
Local object
Stress
Max
Min
Global object
Boundary nodes
図 4.3: 入力を受けて応力が算出された状態
本研究では,物体の破壊は主応力により引き起こされると仮定し,各要素について主応力計算
を行い,全要素の中で主応力が最大となる要素を探索する.最大の主応力を持つ要素の最大主応
力について,破断応力値を超えているか判定を行い,破断応力値を超えていた場合に破壊が進行
する.また,破壊の進行量について記す.全体の主応力に対して閾値を設定し,閾値処理を行う.
37
閾値を超えた領域を破壊進行領域 (Region of fracture progress) として,一度の破壊進行時に,破
壊は破壊進行領域内を直進するものとする.破壊進行領域は,現在のき裂先端から,閾値を超え
た要素の重心位置の中で,最も長い距離を半径とする円状の領域として確保する.図 4.4 に,閾値
処理を行い破壊進行領域が決定された様子を示す.
Region of fracture progress
Stress
Max
Min
図 4.4: 破壊進行領域決定
破断面は,主応力が最大となった要素の最大主応力方向に垂直な面として形成され,破壊が進行
するものと仮定する.破壊が進行した場合,ローカルオブジェクトについてのみ幾何形状と剛性
マトリクスの再構築計算を行う.幾何形状の再構築方法については, 4.2.4 節で詳しく述べる.図
4.5 に,破壊が進行した様子を示す.
Stress
Max
New crack
Min
図 4.5: 破壊進行
破壊進行に伴い,き裂の先端部分を中心としてローカルオブジェクトを更新する.更新する領
域は,現在の破壊進行における破壊進行領域の半径を用いて円状に確保する.図 4.6 に,ローカル
オブジェクトを更新した様子を示す.
38
Stress
Max
Updata of local object
Min
図 4.6: ローカルオブジェクト更新
また,更新の際には,グローバルオブジェクトにローカルオブジェクトのトポロジ変化を反映さ
せるために,ローカルオブジェクトの再構築計算結果をもとに,前処理段階で記録された剛性マ
トリクスの逆行列の中から,類似するデータをグローバルオブジェクトに適用する.類似する逆
行列データというのは,実時間処理中の破壊進行角度と破壊進行量に,最も近い前処理した破壊
進行角度と破壊進行量における逆行列データである.破壊進行回数が i のとき,実時間処理中の破
壊進行角度を θ(i),破壊進行量を L(i) とし,前処理した破壊進行角度を ψj (i),前処理した破壊進
行量を Rk (i), j = 1, 2, · · · , Ncd , k = 1, 2, · · · , Ncl とする. また,実時間処理中の破壊進行角度に
最も近い前処理した破壊進行角度を ψmin (i),破壊進行量を Rmin (i) とする.式 (4.8) に ψmin (i),
式 (4.9) に Rmin (i) を示す.ただし,min は最小値を示す.
ψmin (i) =
Rmin (i) =
min |θ(i) − ψj (i)|
(4.8)
min |L(i) − Rk (i)|
(4.9)
j∈1···Ncd
k∈1···Ncl
類似する前処理した逆行列データは,最初のき裂進展時から現在までの破壊進行角度に最も近い前
処理した破壊進行角度の履歴 ψmin (1), ψmin (2), · · · , ψmin (i) と破壊進行量の履歴 Rmin (1), Rmin (2),
· · · , Rmin (i) によって一意に決定され,物体全体に適用される.
以上の作業を繰り返し行うことで,破壊現象を取り扱う.
39
4.2.3
幾何形状の再構築計算
破壊進行に伴い,幾何構造を再構築する必要がある.本研究では,節点再配置により節点数を
変化させずに幾何形状を変化させる.本節では,幾何形状を取り扱うために必要な基礎的な幾何
学,および簡単な二次元モデルを用いて具体的に幾何形状再構築の過程を示す.
基礎的な幾何学
面 (Plane)S と,二点 A, B を結ぶ線分 AB が交差している場合を考える.面の法線ベクトル
(Normal vector) を n,面上の任意の点の位置ベクトルを X とする.図 4.7 に仮定する状況の概念
図を示す.
Plane S
A
Q
B
X
Normal vector
図 4.7: 交差する面と線分
d を定数,S(t) を線分上の位置として,面は式 (4.10),線分は式 (4.11) で示される.ただし,式
(4.11) 中の t は媒介変数であり,0 ≤ t ≤ 1 が成り立つ.また,OA は原点から点 A へのベクトル,
AB は点 A から点 B までのベクトルである.
d =n·X
(4.10)
S(t) = OA + tAB
(4.11)
面と線分の交点を点 Q としたとき,交点 Q は面 S 上の点であり,線分 AB 上の点でもある.式
(4.11) で示される線分 AB の式 S(t) を,式 (4.10) の X へと代入し,媒介変数 t を求め,得られた
t を式 (4.11) へ代入することで,交点の式が得られる.式 (4.12) に,Q を示す.
d − n · OA
AB
(4.12)
n · AB
また,点 A の面 S への射影距離を k としたとき,式 (4.13) で k が計算される.ただし,AX は
点 A から点 X までのベクトルである.
n (4.13)
k = AX ·
|n| Q = OA +
40
幾何形状再構築手法
簡単な二次元のモデルを用いて,本研究に用いる幾何形状再構築手法 [60] を説明する.本手法
では,物体の節点数を変えることなく,節点間の結合を切り,節点を移動させることで,物体の
トポロジ変化を取り扱う.本手法により,煩雑な四面体の分割 [61]-[65] を行わないために,幾何
形状の再構築計算時間の増加を抑えることが可能である.また,破断面は奥行き方向に一様な平
面であり,一度の進展では,直線的に進展すると仮定する.
破断面と交差する線分を求め,交差線分 (Crossed line) の両端の節点を移動対象節点,交差線分
を削除対象線分,交差線分から構成される要素を削除対象要素とする.図 4.8 に,破断面と交差す
る線分を計算した様子を示す.
Determination of
crossed lines
Crack
Crossed lines
図 4.8: 交差線分算出
次に,移動対象節点を破断面へ射影する.射影距離は,破断面を面 S,削除対象線分を線分 AB
としたとき,式 (4.13) で与えられる.そして,削除対象線分と削除対象要素を取り除くことで,破
断面に平行な形状 (New crack) が形成される.図 4.9 に,節点を射影し,新しい形状が形成された
様子を示す.
Elimination of crossed
lines and crossed elements
New crack
図 4.9: 新たなき裂形状形成の様子
41
4.3
トポロジ変化を詳細に反映した反力提示
前節までに述べた手法では,グローバルオブジェクトについて,破壊によるトポロジ変化を近
似的に反映させるに過ぎず,実時間処理中の物体トポロジを正確に反映した詳細な反力提示は困
難である.詳細な反力提示を行うために,前処理段階で,多様なトポロジと操作についての反力
を記録し,実時間処理中に物体のトポロジとユーザ操作を基に補正して提示する手法を提案する.
本節では,実時間処理段階における反力の補正方法について詳細を記す.実時間処理時には,物
体のトポロジとユーザの操作をもとに,記録された反力データを再生する.しかし,実時間処理時
には,現在のトポロジと,前処理したトポロジが一致しない場合が生じ得るために,提示反力を,
類似するトポロジにおいて記録した反力データで補正して提示する必要がある.以下では,まず,
理想的に現在のトポロジと前処理トポロジが一致する場合について説明する.理想的に現在のト
ポロジと前処理トポロジが一致する場合とは,最初のき裂進展時から現在まで常に前処理された
トポロジを経てきたことを指す.次に,実時間処理中の破壊進展の過程で,一度でも前処理され
ていないトポロジが生じた場合について記す.
4.3.1
理想的に現在のトポロジと前処理トポロジが一致する場合
まず,説明のために記号を導入する.前処理反力データは,破壊進行回数,破壊進行角度,破壊進
行量,ユーザの操作位置,ユーザの操作量,および,ユーザの操作量という六つのパラメータで決
定される.破壊進行回数 i のとき,実時間処理中の破壊進行角度を θ(i),前処理した全 Ncd パター
ンの破壊進行角度を ψ1 (i), ψ2 (i), · · · , ψNcd (i) とする.また,実時間処理中の破壊進行量を L(i),
前処理した全 Ncl パターンの破壊進行量を R1 (i), R2 (i), · · · , RNcl (i) とする.ユーザの操作方向に
ついては,実時間処理中のユーザ操作方向を η ,前処理されたユーザ操作方向を ζ1 , ζ2 , · · · , ζNmd
とする.また,実時間処理中の操作量ベクトルを M s ,前処理したユーザ操作量ベクトルを M p ,
ユーザの操作位置を Pm とする.
理想的に現在のトポロジと前処理トポロジが一致する場合では,破壊進行回数 n のとき,お
よび,過去 n − 1 回を含めた破壊進行角度 θ(1), θ(2), · · · , θ(n) が全て,前処理した破壊進行角度
ψ1 (i), ψ2 (i), · · · , ψNcd (i) に含まれる.さらに,破壊進行量についても同様に,過去 n − 1 回も含
めた破壊進行量 L(1), L(2), · · · , L(n) が全て,前処理した破壊進行量 R1 (i), R2 (i), · · · , RNcd (i) に
含まれる.ただし,i = 1, 2, · · · , n である.仮に,破壊進行回数 i 回目の破壊進行角度 θ(i) が,
前処理した破壊進行角度 ψj (i) と等しいとする.また,破壊進行量 L(i) が,前処理した破壊進
行量 Rk (i) と等しいとする.ただし, i = 1, 2, · · · , n,j = 1, 2, · · · , Ncd ,k = 1, 2, · · · , Ncl で
あるとする.また,実時間処理時の操作方向 η と,前処理した操作方向 ζ1 , ζ2 , · · · , ζNmd の間に,
ζl ≤ η < ζl+1 の関係があるとする.ただし,l = 1, 2, · · · , Nmd である.破壊進行回数 n のときの
提示反力 f n,η,Pm (θ(1), · · · , θ(n), L(1), · · · , L(n)) は,vl , vl+1 を式 (4.14) で表したとき,式 (4.15)
で表される.
vl =
vl+1 =
ζl+1 − η
ζl+1 − ζl
η − ζl
ζl+1 − ζl
42
(4.14)
f n,η,Pm (θ(1), · · · , θ(n), L(1), · · · , L(n)) =
l+1
a=l
va
M sT
M pT
× f n,ζa ,Pm (ψj (1), · · · , ψj (n), Rk (1), · · · , Rk (n))
(4.15)
理想的に現在のトポロジと前処理トポロジが一致する場合では,実時間処理中のトポロジと同
じトポロジについて反力データが記録されているために,記録された反力データに操作方向につ
いての補正を行うだけで反力提示が可能である.
4.3.2
現在のトポロジと前処理トポロジが異なる場合
一回目のトポロジ変化時における反力の補正と,二回目のトポロジ変化時における反力の補正
について記した後,任意の破壊進行回数の場合について記す.
まず,一回目の破壊進行について, 現在の破壊進行角度と前処理した破壊進行角度が一致しな
い状況を考える.まず,実時間処理中のユーザ操作方向 η と,前処理した操作方向 ζ1 , ζ2 , · · · , ζNmd
の間に,ζl ≤ η < ζl+1 の関係があるとする.ただし,l = 1, 2, · · · , Nmd − 1 とする.また,実時
間処理中の破壊進行角度を θ(1) とし,前処理した破壊進行角度とは ψi (1) ≤ θ(1) < ψi+1 (1) の関
係にあるとする.ただし,i = 1, 2, · · · , Ncd − 1 である.さらに,実時間処理中の破壊進行量を
L(1) とし,前処理した破壊進行量とは Rm (1) ≤ L(1) < Rm+1 (1) の関係にあるとする.ただし,
m = 1, 2, · · · , Ncl − 1 である.図 4.10 に,一回目の破壊進行における,現在の破壊進行角度と,前
処理した破壊進行角度の関係を示す.ただし,簡単のため,破壊進行量が同じ場合について示す.
y
Pre-calculated crack direction
Current crack direction
Pre-calculated crack direction
x
図 4.10: 一回目の現在の破壊進行角度と前処理した破壊進行角度の関係
43
本研究では,現在の破壊パターンにおける反力を,類似する前処理破壊進行角度 ψi (1), ψi+1 (1),
類似する破壊進行量 Rm (1), Rm+1 (1),および,類似する操作パターン ζl , ζl+1 における反力の重み
付き和で補正する.破壊進行角度に関する重みを w1,i , w1,i+1 とすると式 (4.16),破壊進行量に関す
る重みを u1,m , u1,m+1 とすると式 (4.17),操作パターンに関する重みを vl , vl+1 とすると式 (4.18)
で示される.
w1,i =
w1,i+1 =
u1,m =
u1,m+1 =
ψi+1 (1) − θ(1)
ψi+1 (1) − ψi (1)
θ(1) − ψi (1)
ψi+1 (1) − ψi (1)
Rm+1 (1) − L(1)
Rm+1 (1) − Rm (1)
L(1) − Rm (1)
Rm+1 (1) − Rm (1)
vl =
vl+1 =
ζl+1 − η
ζl+1 − ζl
η − ζl
ζl+1 − ζl
(4.16)
(4.17)
(4.18)
式 (4.16),式 (4.17),および,式 (4.18) で示した重みを,一回目の破壊進行時における反力補正
重みとして記録し,二回目以降の反力補正にも用いる.一回目の破壊進行時の提示反力は式 (4.19)
で示される.
l+1
i+1 m+1
M sT f 1,η,Pm (θ(1), L(1)) =
vc
w1,a u1,b f 1,ζc ,Pm (ψa (1), Rb (1))
M p T c=l a=i b=m
(4.19)
次に,二回目の破壊進行について考える.実時間処理中の二回目の破壊進行角度を θ(2),前処理
した破壊進行角度との間に ψj (2) ≤ θ(2) < ψj+1 (2) の関係があるとする.また,破壊進行量につい
ても,二回目の破壊進行量を L(2) とし,前処理した破壊進行量との間に Rn (2) ≤ L(2) < Rn+1 (2)
の関係があるとする.ただし,j = 1, 2, · · · , Ncd − 1, n = 1, 2, · · · , Ncl − 1 である.図 4.11 に,二
回目の現在の破壊進行角度と前処理した破壊進行角度の関係を図示する.ただし,二回の破壊進
行において破壊進行量が同じである場合について示す.
44
y
Current crack direction
Pre-calculated crack direction
Past crack direction
x
図 4.11: 二回目の現在の破壊進行角度と前処理した破壊進行角度の関係
二回目の破壊進行時の,破壊進行角度における反力補正重み w2,j , w2,j+1 は,式 (4.20) で示さ
れる.
w2,j =
w2,j+1 =
ψj+1 (2) − θ(2)
ψj+1 (2) − ψj (2)
θ(2) − ψj (2)
ψj+1 (2) − ψj (2)
(4.20)
また,二回目の破壊進行時の,破壊進行量における反力補正重み u2,n , u2,n+1 は,式 (4.21) で示
される.
u2,n =
u2,n+1 =
Rn+1 (2) − L(2)
Rn+1 (2) − Rn (2)
L(2) − Rn (2)
Rn+1 (2) − Rn (2)
(4.21)
提示反力は,一回目の破壊進行角度と反力補正重みも考慮して,式 (4.22) となる.
f 2,η,Pm (θ(1), θ(2), L(1), L(2)) =
l+1
i+1 m+1 j+1 n+1
M sT v
w1,a u1,b w2,c u2,d
e
M p T e=l a=i b=m c=j d=n
× f 2,ζe ,Pm (ψa (1), ψc (2), Rb (1), Rd (2))
(4.22)
最後に,任意の破壊進行回数の場合について記す.現在の破壊進行回数を n として,過去の n − 1
回の破壊進行も含めた破壊進行角度を θ(1), θ(2), · · · , θ(n) とする.また,前処理された破壊進行角
度との間に式 (4.23) の関係が成り立つとする.ただし,i = 1, 2, · · · , n とし,s(i) = 1, 2, · · · , Ncd −1
とする.
ψs(i) (i) ≤ θ(i) < ψs(i)+1 (i)
45
(4.23)
また,過去の n − 1 回の破壊進行も含めた破壊進行量を L(1), L(2), · · · , L(n) とする.また,前
処理された破壊進行量との間に式 (4.24) の関係が成り立つとする.ただし,t(i) = 1, 2, · · · , Ncl − 1
とする.
Rt(i) (i) ≤ L(i) < Rt(i)+1 (i)
(4.24)
破壊進行回数 i 回目の破壊進行角度に対する反力補正重み wi,s(i) , wi,s(i)+1 は式 (4.25) で与えら
れる.
wi,s(i) =
ψs(i)+1 (i) − θ(i)
ψs(i)+1 (i) − ψs(i) (i)
wi,s(i)+1 =
θ(i) − ψs(i) (i)
ψs(i)+1 (i) − ψs(i) (i)
(4.25)
また,破壊進行量に対する反力補正重み ui,t(i) , ui,t(i)+1 は式 (4.26) で与えられる.
ui,t(i) =
ui,t(i)+1 =
Rt(i)+1 (i) − L(i)
Rt(i)+1 (i) − Rt(i) (i)
L(i) − Rt(i) (i)
Rt(i)+1 (i) − Rt(i) (i)
(4.26)
補正された反力ベクトル f n,η,Pm (θ(1), · · · , θ(n), L(1), · · · , L(n)) は,過去 n − 1 回の破壊進行
の重みも考慮して,式 (4.27) となる.
f n,η,Pm (θ(1), · · · , θ(n), L(1), · · · , L(n)) =
s(1)+1
s(n)+1 t(1)+1
t(n)+1
l+1
M sT v
·
·
·
·
·
·
e
M p T e=l a=s(1)
b=s(n) c=t(1)
d=t(n)
× w1,a · · · wn,b u1,c · · · un,d
× f n,ζe ,Pm (ψa (1), · · · , ψb (n), Rc (1), · · · , Rd (n))
(4.27)
式 (4.27) で示したように,破壊進行回数 n のときの反力は,2n+2 個の記録された反力データを
用いて求めることができる.以上の方法で,本研究では実時間処理中のトポロジ変化を詳細に反
映した反力を提示する.
46
第 5 章 システム実装
提案手法に基づき,仮想物体の破壊現象を実時間で取り扱うことができるシステムを構築した.
本章では,まず,システム環境の概要について記す.次に,実装に用いたモデルやパラメータの仕
様について記す.さらに,提案手法の定量的な評価として,従来法との再構築計算時間の比較を
行った結果を記す.最後に,提案手法の反力補正の評価として,反力誤差を調査した結果を記す.
5.1
システム概要
実装システムの全体図を図 5.1 に示す.本システムは,データの入出力を行う力覚提示デバイス
(Novint 社製 Falcon[44]) と表示装置(LCD), および計算機(PC)から構成される.
また,各機器の仕様は表 5.1,実装環境は表 5.2 に示す.
Haptic device
LCD
図 5.1: システム全体図
47
PC
PC
LCD
力覚提示デバイス
表 5.1: 機器の仕様
CPU
Intel 社製 Core 2 Quad 2.4GHz (4 core)
コア総数
4
Memory
4.0GB
GPU
NVIDIA 社製 GeForce8800GTX
製品名
LG 社製 FLATRON L194WT-BF
解像度
WXGA+ 1440×900
製品名
Novint 社製 Falcon
作業空間
101.6mm × 101.6mm × 101.6mm
位置分解能
0.06mm
インタフェース
USB2.0
OS
開発環境
数値計算ライブラリ
メッシュ生成
5.2
5.2.1
表 5.2: 実装環境
Microsoft 社製 Windows XP Professional
Microsoft 社製 Visual Studio 2005
Intel 社製 Math Kernel Library 10.0
Mercury 社製 Amira 4.1
提案モデル実装
モデルの仕様
実装に用いたモデルは,等方性線形弾性体の三次元板形状モデルである.初期状態で,破壊さ
れた領域を有している.図 5.2 に,モデル外観を示す.また,表 5.3 に,実装に用いたモデルの仕
様を示す.
Initial crack
Y Initial local object
X
Z
図 5.2: モデル外観
48
表 5.3: 実装に用いたモデルの仕様
節点数
要素数
大きさ
ヤング率
ポアソン比
破断応力値
破壊進行閾値
5.2.2
506
1266
240mm×240mm×10mm
1.0MPa
0.4
0.4MPa
0.2MPa
前処理パラメータ選定
本研究では,前処理段階で,多様なトポロジにおける剛性マトリクスの逆行列と,多様なトポロ
ジとユーザ操作パターンにおける反力を記録する.しかし,全てのトポロジやユーザ操作に対し
て記録することは不可能であるために,パターンを限定して記録する.前処理パラメータは,前
章でも示したように,下記のようになる.
• 破壊進展数 (Ncp )
• 破壊進行角度 (Ncd )
• 破壊進行量 (Ncl )
• ユーザの操作点 (Nmp )
• ユーザの操作方向 (Nmd )
• ユーザの操作量 (Nml )
上記パラメータの中で,破壊進展数と破壊進行角度,および,破壊進行量は,剛性マトリクス
の逆行列と反力で共通のパラメータであり,ユーザの操作点,操作方向,操作量は反力限定のパラ
メータである.まず,共通パラメータである破壊進展量と破壊進行角度,破壊進行量に関して記
す.破壊進展数について,前章でも示したように前処理で記録するデータ量は,破壊進展数乗の
オーダで増大する.また,本手法では,実時間処理中に前処理した剛性マトリクスの逆行列デー
タは,実メモリ上に配置することからも,残りのパラメータを全て決定した後に,計算機のメモ
リ容量から逆算して決定する.また,破壊進行角度が異なるモデルは剛性マトリクスが異なるた
めに,同じ操作をした場合の,変形や反力も異なる.破壊進行角度が異なる二つのモデルを考え
たとき,破壊進行角度が大きく異なる場合は,同じ操作に対する変形や反力が大きく異なると考
えられる.また,逆に,破壊進行角度があまり変わらない場合は,同じ操作をしたときの変形や反
力の違いは小さいと考えられる.本研究では,ヒトは,どの程度の破壊進行角度の違いを力覚的
に分別可能であるかを被験者実験によって調査し,破壊進行角度を選定する目安とした.被験者
実験については,5.2.3 で詳しく記す.破壊進行量については,破壊進行角度と同様に,ヒトの力
覚に大きく影響を与えるパラメータである.ただし,本研究で扱う線形有限要素モデルでは,要
素内の応力は一定となるために,破壊進行量については一要素以上の分解能が得られない.一方,
破壊進行角度は,応力値に影響されるために,要素内の応力が一定であるとしても,値が異なれ
49
ば多様なパターンが発生する.従って,本研究では,より多様なパターンを発生させるパラメー
タである破壊進行角度数を確保することを優先し,破壊進行量は一通りとした (Ncl = 1).なお,
破壊進行量として,初期状態のモデルの全自由節点に対して,後述するユーザ操作量と操作方向
の全ての組み合わせで,変位を与えたときに算出される破壊進行量の平均値を用いた (40.3mm).
ユーザ操作量のパターン数については,線形モデルでは変位と外力間に線形関係があると仮定
できるため,1 パターンのみとした (Nml = 1).操作量は,モデルの奥行き方向に垂直な面の一辺
の長さの 12 である 120mm とした.ユーザの操作点数は,ユーザが引っ張る点は自由表面節点に限
られるとして,全自由表面節点とした.表 5.3 に示したモデルの場合,自由表面節点の総数は 428
である (Nmp = 428).ユーザの操作は,板状のモデルを引っ張り破壊する操作を仮定し, xy 平面
内に限定した.また,操作方向は,面内操作で可動な 360 度の範囲を,10 度おきに 36 分割し,36
パターンとした (Nmd = 36).
5.2.3
被験者実験
前処理する破壊進行角度を選定する目安として,人間の破壊進行角度に対する力覚的な分解能
を被験者実験により調査した.破壊進行角度に対する力覚的な分解能とは,力覚提示デバイスを
用いて破壊進行角度の異なるモデルを操作したときに,どの程度の破壊進行角度の違いであれば
力覚的に判別できるかという閾値である.図 5.3 に,被験者実験の様子を示す.
Examinee
Display
Haptic device
図 5.3: 被験者実験の様子
被験者実験では,各試行において見本とモデル A,B という三つのモデルを提示した.図 5.4 に,
見本とモデル A,B の例を示す.見本は,初期き裂から破壊が x 軸から 90deg の角度へ進行したモ
デルであり,モデル A と B のうち一方は,見本と同一のモデル,他方は見本と異なる角度 θcurrent
に破壊が進行したモデルである.ただし,0deg< θcurrent < 90deg である.
50
Sample
Model A
90deg
Model B
90deg
(a) (b) (c)
図 5.4: 被験者実験に用いた見本とモデル A,B
(a) 見本 (b) モデル A (c) モデル B
被験者には,モデル A と B のどちらが見本と同じモデルであるか伝えずに提示した.提示順序
は,最初に見本を提示し,モデル A と B の提示順序は無作為とした.操作は,力覚提示デバイス
を用いて行い,仮想空間内のモデルを把持し,引っ張る操作である.被験者には,視覚的には常
に見本について,入力に応じた変形計算の結果を提示し,力覚的にはモデル A,および,モデル
B を操作した時の反力を提示した.視覚情報を提示したことにより,被験者は視覚的にも自身の操
作点移動量を確認することができる.ただし,実験条件を統一するために,接触点は二通り,操
作方向は三通りに限定した.接触点は,仮想空間にどこに操作点があっても,指定した点のみを
把持することで実現した.また,操作方向は,デバイスからの入力の中で,指定した方向成分の
みを処理することで実現した.実際には,実験時に被験者に矢印で操作方向を示している.また,
作用させる変位の大きさは,被験者の自由として制限せずに実験を行った.
見本とモデル A,B に対する操作が終わったあと,被験者にモデル A と B のうち,どちらが見
本と同じモデルであると感じたかを質問した.正解の場合,次に比較する方向 θnext を式 (5.1) の
ように,90deg へ近づけた.
θnext = θcurrent +
(90 − θcurrent)
2
(5.1)
不正解の場合は,次に比較する方向 θnext を式 (5.2) のように,90deg から遠ざけた.
θnext = θcurrent −
(90 − θcurrent)
2
(5.2)
ただし,回答の正解・不正解は被験者に伝えないこととした.上記の作業を初期進行角度を 45deg
として 5 回繰り返して行い,操作点と操作方向の組み合わせ毎に 2 回ずつ行った.被験者は,右
利きの男性 7 名である.
表 5.4 に,破壊進展方向毎の比較回数と正解数,および,1:2 点比較法に基づいた検定におけ
る有意水準 5 パーセントでの棄却限界を示す.棄却限界とは,結果が有意かどうか判別するため
の限界値であり,正解数が棄却限界を上回った場合は,有意に判別可能と言える.逆に,正解数
が棄却限界を下回った場合は,有意に判別可能であるとは言えない.また,棄却限界値は試行回
数毎に異なる値を持つために,表 5.4 では,各度数の試行回数毎に棄却限界値を示しており,各度
数の正解数が棄却限界を超えた場合には,正解数の右上に*を付して示している.
51
表 5.4: 被験者実験の結果
破壊進行角度の差 (deg)
試行回数
正解数
棄却限界
0 度以上 15 度未満
15 度以上 30 度未満
30 度以上 45 度未満
45 度以上 60 度未満
60 度以上
124
128
41
114
14
58
85*
29*
92*
13*
72
74
27
67
11
実験の結果,破壊進行角度の差が 15deg 以上の場合は,いずれの場合も棄却限界以上の正解数
が得られ,破壊進行角度の差が有意に違いが判別できていると言える.しかし,0deg 以上 15deg
未満の場合は,正解数が棄却限界を下回り,違いを判別できているとは言えなかった.
したがって,本実験により区別できる最小の破壊進行角度の差は,15deg とした.なお,過去の
研究や様々な文献において,ヒトの破壊進行角度に対する力覚的分解能についての調査は見当た
らず,本実験結果は,本実験で用いたモデル限定ではあるものの,ヒトの知覚特性について新し
い知見を与える結果であると言える.
前処理パラメータの設定では,前処理した破壊進行角度と実時間処理時の破壊進行角度の差が,
15deg 以内となるように設定した.また,破壊は前に進むという仮定をおき,破壊進行角度がとり
得る範囲を 0deg から 180deg までに限定した.結果,前処理する破壊進行角度は,180deg を 15deg
おきに 13 パターンとして計算した.
最後に,破壊進行回数を決定した.前処理データの中で,剛性マトリクスの逆行列はメモリ上に
配置し,反力データはハードディスクに記録する.本研究で用いた計算機のメモリ容量は 4GByte
であり,前述したパラメータについて全てのパターンの逆行列データを記録するとして考えると,
記録できる最大の破壊進行回数を x として,剛性マトリクスの逆行列一つあたりの容量 18MByte
を考慮すると,前章の式 (4.1) を用いて,式 (5.3) を満たすように,x を選ぶ必要がある.
6
(18.0 × 10 × 1) ×
x
13i ≤ 4.0 × 109
(5.3)
i=1
式 (5.3) から,x の値を求めると,x ≤ 2 となり,記録できる最大の破壊進展回数は 2 回である.
表 5.5 に,本研究で採用した前処理パラメータを示す.
表 5.5: 前処理パラメータ
破壊進行回数
2回
13 通り
破壊進行量
1 通り (40.3mm)
ユーザ操作点
428 通り
ユーザ操作方向
36 通り
ユーザ操作量
1 通り (120.0mm)
データ量(逆行列)
約 3.3GByte
データ量(反力)
約 67.3MByte
破壊進行角度
52
また,前処理では大量の逆行列計算や,反力計算を行うために,汎用の PC では過度の負担がか
かる.実装にあたっては,大規模な計算を想定して設計されているクラスタ計算機を用いて行っ
た.表 5.6 に,クラスタ計算機の仕様,図 5.5 にクラスタの外観を示す.
クラスタ計算機
LCD
表 5.6: クラスタ計算機の仕様
CPU
AMD 社製 Opteron 2.0GHz (4 core) ×4
コア総数
16
Memory
32GB
GPU
NVIDIA 社製 Quadro FX 5500
製品名
BenQ 社製 FP2092
解像度
UXGA 1600×1200
Display Haptic device
Cluster
Cluster
(a) (b) 図 5.5: クラスタ外観
(a) 全体図 (b) クラスタ拡大図 53
5.3
破壊操作の結果例
本手法を用いて,ユーザの操作によって仮想物体を破壊している様子を図 5.6 に示す.図 5.6 で
は,ローカルオブジェクトは白いポリゴン,グローバルオブジェクトは灰色のポリゴンで示して
いる.
Local object
Reaction force vector
図 5.6: ユーザ操作による仮想物体の破壊結果
提案手法により,ユーザの操作によって仮想物体を破壊可能であることが示された.本手法を
用いることにより,仮想物体の弾性変形のみならず,従来扱うことが困難であった破壊現象を扱
うことができる.
54
5.4
定量評価
本手法の定量的な評価として,従来法との再構築時間の比較,および,提示反力の比較を行った.
5.4.1
再構築計算時間比較
破壊進行時に,物体全体の剛性マトリクスの逆行列の再構築計算を行った場合の計算時間と,提
案手法により,ローカルオブジェクトのみの再構築計算を行った場合の計算時間を比較した.計
測は表 5.1 に示した機器で行った.入力として,全てのモデルに同一の変位を与え,破壊進行時に
おける再構築時間の平均値を,各モデルにおける再構築時間とした.比較に用いたモデルの仕様
を表 5.7 に,再構築計算時間測定結果を図 5.7 に示す.なお,全モデルについて,一定の変位を与
えたとき,2 回の破壊進行が生じた.
表 5.7: 再構築時間比較に用いたモデルの仕様
節点数
要素数
506
1266
765
2017
1012
2718
1270
3487
6000
Proposed method
5905.7
Calculation time (ms)
5000
Conventional method
4000
2813.1
3000
2000
1290.9
376.9
1000
12.6
7.4
30.4
27.5
0
0
500
750
1000
1250
The number of nodes
図 5.7: 計算時間測定結果
提案手法により,従来法に比べて再構築計算時間が削減された.また,節点数が多くなるほど,
計算時間削減の効果が大きくなることが示された.
55
5.4.2
提示反力の評価
4.3 節で記したように,提案手法では,前処理した反力データを補正してユーザに提示する.反
力データを補正せずに提示した場合と,提案手法のように補正した場合について,ローカルオブ
ジェクトを用いない従来法における反力を正しい値として,誤差を比較した.比較に用いたパラ
メータを表 5.8 に示す.
表 5.8: 提示反力誤差を調査する際に用いたパラメータ
進行回数
進行角度
操作点
操作方向
2回
180 通り
428 点
36 方向
各操作点について,操作方向毎に全進行パターンにおける反力誤差ベクトルを,反力補正あり
の場合と補正無しの場合について求め,反力誤差ベクトルの大きさを,該当操作点における反力
誤差とした.また,反力誤差の程度を評価するために,反力誤差の平均値と標準偏差を算出した.
ただし,破壊進行回数毎に,反力誤差を正規化した.進展回数 1 回における結果を表 5.9,進展回
数 2 回における結果を表 5.10 に示す.
表 5.9: 提示反力評価結果(進展回数 1 回)
2773440
2574967
92.8 %
−3
4.0 × 10 ± 2.6 × 10−2
2.1 × 10−3 ± 1.9 × 10−2
評価パターン数
誤差が小さくなった数
誤差が小さくなった割合
平均誤差 (補正無し)
平均誤差 (補正あり)
表 5.10: 提示反力評価結果(進展回数 2 回)
499219200
355727627
71.3 %
−3
4.8 × 10 ± 2.7 × 10−2
3.2 × 10−3 ± 1.7 × 10−2
評価パターン数
誤差が小さくなった数
誤差が小さくなった割合
平均誤差 (補正無し)
平均誤差 (補正あり)
評価を行った多くのパターンについて,提案手法における反力補正により,補正を行った場合
の方が,補正を行わない場合よりも,真値との誤差が小さくなったことが示された.また,誤差
の程度を調査した結果,いずれの進展回数においても,補正を行った場合の方が,補正を行わな
い場合よりも,誤差の平均値と標準偏差が小さくなった.提案手法における反力補正は有効であ
ると言える.
ただし,一部のパターンについては,反力補正を行った場合の方が,補正を行わない場合より
56
も誤差が大きくなることも確認された.また,平均値や標準偏差についても,補正を行った場合の
方が,補正を行わない場合に比べて低いというものの,オーダとしては同じであった.提案した
補正方法は,大半の操作やトポロジには有効であると言えるものの,不十分な手法であることも
事実である.誤差が増加した要因としては,破壊進行角度の違いによる反力の違いに,線形関係
を仮定して補正したことが考えられる.提案手法では,実時間処理時の破壊進行角度が線形的に
増加するにつれて,前処理反力データに乗算する重みは線形的に増加,もしくは,減少する.提
案手法では,破壊進行角度が線形的に増加した場合に,反力値が線形的に変化すると仮定してい
るためである.しかし,評価において,一部で補正を行った場合の方が,補正を行わない場合に
比べて誤差が大きくなるパターンが存在することから,破壊進行角度の変化と反力の変化に線形
関係がないパターンが存在すると考えられる.
57
第 6 章 おわりに
本研究では,従来では困難であった仮想物体の破壊現象を,実時間で扱う手法を提案した.提
案手法は,再構築時間の削減手法,および,破壊進行時の反力提示手法で構成される.再構築時
間の削減は,毎回の破壊進行時に,物体の破壊領域近傍の局所領域についてのみ再構築計算のみ
を行うことで実現した.また,実時間処理中のトポロジ変化を詳細に反映した反力を提示するた
めに,反力の前処理段階での記録,および,実時間処理中の補正を行った.
前処理する際のパラメータは,被験者実験によってヒトの破壊進行角度に対する力覚的分解能
を調査し,破壊進行角度パラメータの場合は,最大の角度差が最大で 15deg となるように選定す
ればよいと判明した.また,提案手法による前処理と実時間処理を行うことで,従来の有限要素
モデルでは扱うことが困難であった破壊現象を実時間で扱うことができた.また,幾何形状と剛
性マトリクスの再構築計算時間を従来法と比較した結果,提案手法により計算時間が削減された.
さらに,提案手法による反力補正の有効性を示すために,補正を行う場合と補正を行わない場合
の誤差を比較した.結果,評価した多くのパターンにおいて,補正を行った場合の方が,補正を
行わなかった場合に比べて,従来法における反力との誤差が小さいことがわかり,提案手法によ
る補正が有効であることが示された.
本研究では,簡単のため,単純なモデルを用いたが,提案した手法自体は,任意の三次元モデ
ルについても適用可能である.今後は,本手法が適用された力覚提示 VR システムが,破壊現象
も含む操作を扱うことが可能となることで,弾性変形だけでは不可能であった様々な分野・用途
に応用されることが期待される.
58
謝辞
本研究は,大阪大学大学院基礎工学研究科 大城研究室で行ったものである.本研究を遂行する
にあたり,研究の機会を与えてくださった大阪大学大学院基礎工学研究科 大城 理 教授に篤く御
礼申し上げます.修士論文の執筆にあたり,副査として数多くの貴重,かつ,有益な御助言を頂
きました大阪大学大学院基礎工学研究科 河原 源太 教授に篤く御礼申し上げます.研究内容につ
いて多くの御指摘を頂きました大阪大学大学院基礎工学研究科 黒田 知宏 准教授に感謝致します.
研究方針に関して多くの助言を頂きました大阪大学大学院基礎工学研究科 黒田 嘉宏 助教に感謝
致します.修士二年時から,丁寧に研究指導を頂きました大阪大学臨床医工学融合研究教育セン
タ 鍵山 善之特任助教に感謝致します.さらに,修士過程の二年間,共に助け合い励まし合って
きた大阪大学大学院基礎工学研究科 大城研究室の皆様に御礼申し上げます.なお,本研究におけ
る被験者実験は,大阪大学大学院基礎工学研究科人を対象とした研究に関する倫理委員会の承認
(20-6) を得て行った.
59
業績
・国内学会
荒井 良祐,黒田 嘉宏,黒田 知宏, 大城 理. ” 仮想生体膜のき裂進展シミュレーション”,平成
19 年度 電気関係学会関西支部連合大会,神戸 (2007).
荒井 良祐,黒田 嘉宏,黒田 知宏, 大城 理. ” 破壊力学的手法を用いた生体軟組織の破壊進行シ
ミュレーション”,第 47 回日本生体医工学会大会,神戸 (2008).
荒井 良祐,黒田 嘉宏,鍵山 善之,黒田 知宏, 大城 理. ” 応力状態による破壊進行を考慮した
有限要素モデルの構築”,第 13 回日本バーチャルリアリティ学会大会,生駒 (2008).
荒井 良祐,黒田 嘉宏,鍵山 善之, 黒田 知宏, 大城 理. ” 応力状態を考慮した破壊操作が可能な
有限要素モデルの構築”, ME とバイオサイバネティックス研究会,四条畷 (2008).
荒井 良祐,黒田 嘉宏,鍵山 善之, 黒田 知宏, 大城 理. ” トポロジ変化の実時間力視覚レンダリ
ング”,第 53 回システム制御情報学会研究発表講演会,神戸 (2009), 投稿中.
・その他
荒井 良祐,黒田 嘉宏,黒田 知宏, 大城 理. ” 仮想生体膜のき裂進展シミュレーション”,ジョ
イントセミナ,豊中 (2007).
荒井 良祐,黒田 嘉宏,黒田 知宏, 大城 理. ” 破壊進行を考慮した有限要素モデルの構築”,河
原研・大城研 合同ゼミ,豊中 (2008).
・受賞
平成 19 年度 電気関係学会関西支部連合大会奨励賞
60
参考文献
[1] N.Molino, Z.Bao and R.Fedkiw, ”A Virtual Node Algorithm for Changing Mesh Topology
During Simulation”, ACM Transactions on Graphics, Vol.23, No.3, pp.385-392, 2004.
[2] U.Kuhnapfel, H.K.Cakmak, H.Maab, ”Endoscopic Surgery Training Using Virtual Reality
and Deformable Tissue Simulation”, Computer and Graphics, Vol.24, No.5, pp.671-682,
2000.
[3] J.Berkley, G.Turkiyyah, D.Berg, M.Ganter and S.Weghost,”Real-time Finite Element Modeling for Surgery Simulation: an Application to Virtual Suturing”,IEEE Transactions on
Visualization and Computer Graphics”, Vol.10, No.3, pp.314-325, 2004.
[4] M.Nakao, K.Minato, T.Kuroda, M.Komori, H.Oyama, T.Takahashi, ”Transferring Bioelasticity Knowledge through Haptic Interaction”, IEEE MultiMedia, Vol.13, No.3, pp.50-60,
2006.
[5] J.Berkley, S.Wedghorst, H.Gladstone, G.Raugl, D.Berg, M.Ganter, ”Banded Matrix Approach to Finite Element Modeling for Soft Tissue Simulation”, Virtual RealityResearch
Development and Application, Vol.4, No.3, pp.203-212, 1999.
[6] N.Mukai, M.Harada, K.Muroi, T.Hikichi and A. Yoshida, ”New Graphics Models for PC
Based Ocular Surgery Simulator”, Proc of Medicine Meets Virtual Reality 2001, IOS Press,
pp.329-335, 2001.
[7] M.Sedef, E.Samur and C.Basdogan, ”Real-Time Finite-Element Simulation of Linear Viscoelastic Tissue Behavior Based on Experimental Data”, IEEE Computer Graphics and
Applications, Vol.26, No.6, pp.58-68, 2006.
[8] C.Basdogan, M.Sedef, M.Harders and S.Wesarg, ”VR-Based Simulators for Training in
Minimally Invasive Surgery”, IEEE Computer Graphics and Applications, Vol.27, No.2,
pp.54-66, 2007.
[9] C.Basdogan, S.De, J.Kim, M.Muniyandi, H.Kim and M.A.Srinivasan, ”Haptics in Minimally Invasive Surgical Simulation and Training”, IEEE Computer Graphics and Applications, Vol.24, No.2, pp.56-64, 2004.
[10] 黒田 嘉宏,平井 真,中尾 恵,佐藤 寿彦,黒田 知宏,長瀬 啓介,吉原 博幸,” 多指力覚提
示装置を用いた臓器圧排シミュレータに関する研究” 日本バーチャルリアリティ学会論文誌,
Vol.11,No.4,December 2006,pp.515-525
[11] K.Hirota and T.Kaneko, ”Haptic Representation of Elastic Objects”, Presence: Teleoperators and Virtual Environments, Vol.10, No.5, pp.525-536, 2001.
61
[12] S.P.DiMaio and S.E.Salcudean, ”Interactive Simulation of Needle Insertion Models”, IEEE
Transactions on Biomedical Engineering, Vol.52, No.7, pp.1167-1179, 2005.
[13] S.P.DiMaio and S.E.Salcudean, ”Needle Insertion Modeling and Simulation”,IEEE Transactions on Robotics and Automation, Vol.19, No.5, pp.864-875, 2003.
[14] S.P.DiMaio and S.E.Salcudean, ”Needle Steering and Motion Planning in Soft Tissues”,
IEEE Transactions on Biomedical Engineering, Vol.52, No.6, pp.965-974, 2005.
[15] 本間 達,若松 秀俊,” ノコギリによる切離抗力の表現”, 電子情報通信学会論文誌, Vol.J84-A,
No.6, pp.860-869, 2001.
[16] 広田 光一,田中 厚子,金子 豊久,” 力覚をともなう変形・切断操作”, 日本機械学会ロボティ
クス・メカトロニクス講演会’99 講演論文集, 2P1-34-046(1)-2P1-34-046(2), 1999.
[17] W.Wu and P.A.Heng, ”A hybrid condenced finite element model with GPU acceleration
for interactive 3D soft tissue cutting”, Computer Animation and Virtual Worlds, Vol.15,
No.3-4, 2005.
[18] M.Nakao, T.Kuroda, H.Oyama, G.Sakaguchi and M.Komeda, ”Physics-Based Simulation
of Surgical Fields for Preoperative Strategic Planning”, Journal of Medical Systems, Vol.30,
No.5, pp.371-380, Octorber 2006.
[19] D.Bielser, P.Glardon, M.Teschner, M.Gross, ”A state machine for real-time cutting of tetrahedral meshes”, Graphical Models, Vol.66, No.6, pp.398-417, 2004.
[20] A.Lindblad and G.Turkiyyah, ”A Physically-Based Framework for Real-time Haptic Cutting and Interaction with 3D Continuum Models”, Proceedings of the 2007 ACM symposium
on Solid and physical modeling, pp.421-429, 2007.
[21] D.Bielser, V.A.Maiwald and M.H.Gross, ”Interactive Cuts through 3-Dimensional Soft Tissue”, EUROGRAPHICS ’99, Vol.18, No.3, pp.31-38, 1999.
[22] D.Steinemann, M.Harders, M.Gross and G.Szekely, ”Hybrid Cutting of Deformable Solids”,
Virtual Reality Conference, pp.25-29, 2006.
[23] D.Bielser and M.H.Gross, ”Interactive Simulation of Surgical Cuts”, Proceedings of Pacific
Graphics 2000, pp.116-125, 2000.
[24] W.Wu and P.A.Heng, ”An Improved scheme of an interactive finite element model for 3D
soft-tissue cutting and deformation”, Visual Computer, Vol.21, No.8-10, pp.707-716, 2005.
[25] M.Mahvash, ”Novel Approach for Modeling Separation Forces Between Deformable Bodies”, IEEE Transactions on Information Technology in Biomedicine, Vol.10, No.3, 2006.
[26] D.Kobayashi, S.Mizuno, J.Toriwaki and S.Yamamoto, ”Virtual Sculpting with a Pressure
Sensitive Pen”,Proceeding of ACM SIGGRAPH 2003 Sketch & Applications, 2003.
[27] H.Chen and H.Sun, ”Real-time Haptic Sculpting in Virtual Volume Space”, Proceeding of
the ACM symposium on Virtual reality software and technology, pp.81-88, 2002.
62
[28] B.Baxter, V.Scheib, C.Lin and D.Manocha, ”DAB: Interactive Haptic Painting with 3D
Virtual Brushes”,Proceeding of ACM SIGGRAPH 2001, pp.461-468, 2001.
[29] S.Hasegawa, M.Sato, Y.Dobashi, T.Yamamoto, M.Kato and T.Nishita, ”Virtual canoe:
real-time realistic water simulation for haptic interaction”, International Conference on
Computer Graphics and Interactive Techniques ACM SIGGRAPH 2005 Emerging technologies, No. 28, 2005.
[30] S.Hasegawa,N. Fujii,M. Sato, ”Real-time Rigid Body Simulation for Haptic Interactions
Based on Contact Volume of Polygonal Objects”, Computer Graphics Forum, Vol.23, No.3,
pp.529-538, 2004.
[31] H. Mitake,S. Hasegawa,Y. Koike,M. Sato, ”Reactive Virtual Human with Bottom-up and
Top-down Visual Attention for Gaze Generation in Realtime Interactions”, IEEE Virtual
Reality 2007, pp.211-214, 2007.
[32] T.Aoki, K.Asano, R.Ayukawa, S.Hasegawa, H.Ichikawa, Y.Iio, T.Kawase, T.Kuriyama,
I.Matsumura, T.Matsushita, H.Mitake, M.Sato, and T.Toyama, ” Kobito -Virtual
Brownies-”, SIGGRAPH 2005 Emerging technologies, Los Angeles USA, July. 2005.
[33] 青木 孝文,三武 裕玄,浅野 一行,栗山 貴嗣,遠山 喬,長谷川 晶一,佐藤 誠,” 実世界で
存在感を持つバーチャルクリーチャの実現 Kobito -Virtual Brownies-”, 日本バーチャルリア
リティ学会論文誌,Vol.11, No.2, pp.313-322, 2006.
[34] 三好 俊郎, ” 有限要素法入門”, 培風館,東京,1994.
[35] L.M.Smith and D.V.Griffiths, ”Programming the Finite Element Method”,JOHN WILEY
& SONS Inc. NewYork, 1998.
[36] 山田 貴博,” 計算力学レクチャーシリーズ9 高性能有限要素法”, 丸善,東京,2007.
[37] 小林 英男,” 破壊力学” 共立出版,東京,1993.
[38] Q.Zhu, Y.Chen and A.Kaufman, ”Real-time Biomechanically-based Muscle Volume Deformation using FEM”, Computer Graphics Forum, Vol.17, No.3, pp.275-284, December
2001.
[39] A.Faraci, F.Bello and A.Darzi, ”Soft Tissue Deformation using a Nonlinear Hierarchical
Finite Element Model with Real-Time Online Refinement”, Medicine Meets Virtual Reality
13 : The Magical Next Becomes the Medical Now, Vol.111, pp.137-144, 2005.
[40] O.Comas, Z.A.Taylor, J.Allard, S.Ourselin, S.Cotin and J.Passenger, ”Efficient Nonlinear
FEM for Soft Tissue Modelling and Its GPU Implementation within the Open Source
Framework SOFA”, Biomedical Simulation: 4th International Symposium 2008, London,
Uk, pp.28-39, July, 2008.
[41] Z.A.Taylor, M.Cheng and S.Ourselin, ”High-Speed Nonlinear Finite Element Analysis for
Surgical Using Graphics Processing Units”, IEEE Transactions on Medical Imaging, Vol.27,
No.5, May, 2008.
63
[42] J.R.Navarro and A.Susin, ”Non Structured Meshes for Cloth GPU Simulation using FEM”,
Workshop on Virtual Reality Interaction and Physical Simulation, 2006.
[43] SensAble 社 HP:http://www.sensable.com/
[44] Novint 社 HP:http://home.novint.com/
[45] 長谷川 晶一,赤羽 克仁,馬場 次郎,小池 康晴,佐藤 誠,” 高解像度力覚インタフェースを
もつ物理 VR システムの開発”,電子情報通信学会論文誌, J88-D-1(2), pp.1-8, 2005.
[46] 佐藤 誠,” バーチャルリアリティの基礎2/人工現実感の設計:究極のインタフェースを求め
て”,培風館,東京,2000.
[47] Immersion 社:http://www.immersion.com
[48] S.Fukamachi, H.Kajimoto, N.Kawakami and S.Tachi,”Gravity Grabber:Wearable Haptic
Display to present Virtual Mass Sensation”, 34th Int. Conf. On Computer Graphics and
Interactive Techniques(ACM SIGGRAPH 2007), Emerging Technologies, San Diego, USA,
2007.
[49] Y.Tanaka and T.Kanamori,”Dynamic Force Display Device by Pneumatic Pressure Feedback”, Proceeding. FLUCOME, Vol.2, pp.719-723, 1997.
[50] 本間 達,若松 秀俊,” リアルタイムで切離可能な仮想粘弾性物体の構築”, 日本バーチャルリ
アリティ学会論文誌, Vol.6, No.2, pp.137-143, 2001.
[51] 本間 達,若松 秀俊,張 暁林,” 粘弾性体モデルで表現した仮想物体の構築と加工”, 電気学
会論文誌 C,Vol.119, No.12, pp.1437-1443, 1999.
[52] 本間 達,若松 秀俊,” 力覚提示システムを前提とした粘弾性体のナイフによる切離モーメン
トの表現”, 計測自動制御学会論文集, Vol.40, No.4, pp.137-143, 2004.
[53] 本間 達,若松 秀俊,” 力覚提示システムに応用可能なノコギリの切離モーメントの算出”, 日本バーチャルリアリティ学会論文誌, Vol.9, No.3, pp.319-326, 2004.
[54] 本間 達,若松 秀俊,” 粘弾塑性体モデルで表現した物体間の相互作用による破壊”, 計測自動
制御学会論文集, Vol.44, No.7, pp.600-608, 2008.
[55] N.Kume, ”Distributed Massive Simulation for Haptic Virtual Reality Based Surgical Skill
Transfer”, PhD thesis, Kyoto University, March 2006.
[56] 遠藤 恒史,田村 信彦,中口 俊哉,津村 徳道,三宅 洋一, ” 人体腰部の反力測定・再現シ
ステムの構築と腰椎穿刺トレーニングシステムへの応用”, 電子情報通信学会技術研究報告,
Vol.104, No.318, pp.41-45, 2004.
[57] 田川 和義,広田 光一,廣瀬 通孝,” 力覚インタラクションのための動的変形モデル”, 電子
情報通信学会論文誌, Vol.J90-D, No.9, pp.2615-2623, 2007.
[58] K.Tagawa, K.Hirota and M.Hirose,”Impulse Response Deformation Model : an Approach
to Haptic Interaction with Dynamically Deformable Object”, Symposium on Haptic Interfaces for Virtual Environment and Teleoperator Systems 2006,pp.209-215, 2006.
64
[59] S.Saga, K.Vlack, H.Kajimoto and S.Tachi, ”Haptic Video”, 32nd International Conference
On Computer Graphics and Interactive Techniques (SIGGRAPH2005), Los Angeles, USA,
2005.
[60] 中尾 恵,河本 敏孝,杉浦 忠男,湊 小太郎,” 弾性変形モデルに対する頂点数を保存した切
開方法”, 日本バーチャルリアリティ学会論文誌, Vol.12, No.4, pp.585-594, 2007.
[61] P.Alliez, D.Cohen-Steiner, M.Yvinec, and M.Desbrun.”Variational tetrahedral meshing”,
ACM Transactions on Graphics, Vol.24, No.3, pp.617-625, 2005.
[62] F.Ganovelli, P.Cignoni, C.Montani and R.Scopigno, ”Enableing cuts on multiresolution
representation”, The Visual Computer, Vol.17, pp.274-286, 2001.
[63] L.Chen,”Mesh smoothing schemes based on optimal Delaunay triangulations”, In Proceedings of the 13th International Meshing Roundtable, pp.109-120, 2004.
[64] T.Berker,”Automatic mesh generation for complex 3-dimentional regions using a constrained Delaunay triangulation”, Engineering with Computers, Vol.3, No.3, pp.161-175,
1989.
[65] C.Forest, H.Delingette and N.Ayache, ”Removing tetrahedra from manifold tetrahedralisation : application to real-time surgical simulation”, Medical Image Analysis, Vol.9, pp.113122, 2005.
65
ダウンロード

トポロジ変化の実時間レンダリング