リアルタイムシミュレーションとは インターラクティブな作業（ゲーム、ＶＲ）における結果表

示にシミュレーションを用いる
予め用意した結果を表示する→多様な現象を表現できない
 その場でリアルタイムに計算した結果を表示する


CGのためCGのためののためのリアルタイム大変形
リアルタイム大変形シミュレーション
大変形シミュレーション

 計算時間＜実時間が絶対の条件

東京大学 人工物工学研究センター
鈴木克幸

物理ベースＣＧに求められる技術

研究背景 変形シミュレーション
 変形のシミュレーション

 計算時間に対する要求は強い

クリエーターの想像力を妨害しないために、早く結果がほしい できればリアルタイム
 精度は求めない

目に見えない差異はどうでもいい 実際の現象とあっていなくてもいい

バネマスモデル → 大変形までの定式化が簡単
計算量が少ないためリアルタイムで計算可能連続体力学との整合性がとりにくい 有限要素法 →連続体力学により定式化される計算が大規模になりがち

研究背景 リアルタイムシミュレーション
 リアルタイムシミュレーションにおいては、実際にシミュレーションを行なう前にある程度の計算をしておくことが出来る（前処理）
 前処理には現実的な時間であれば時間を掛けてもよい 線形計算では、インパルス応答の線形和などが使われる。

研究背景 変形シミュレーション
線形有限要素法 → 結果を重ね合わせで表現できる
大変形では不自然な結果になる
 非線形有限要素法 → 大変形でも自然な結果になる
重ね合わせによる表示が使えない


MODEL ORDER REDUCTION (MOR)

MODEL ORDER REDUCTION 概要

 近年注目されている手法で、モデルを低次元化するこ Model Order Reduction(以下：MOR)では変位ベクトルをとにより計算時間を高速化する技術
 大規模な設計に対して用いられる高速演算手法であり設計初期段階に解析対象の特性などを短時間で知ることができる手法である
 未知数や自由度の数を大幅に落としながらも、入出力のサイズはそのまま維持できるのが特徴
 制御、流体力学、固体力学など、幅広い分野で研究次の式で近似する。

u = Tξ
= [q1 q 2 ξ1 ξ L qm ]  2 M  ξ m       =  q1  ξ1 +  q 2  ξ 2 + L +  q m  ξ m      uT ξm×1n×1n×m

MODEL ORDER REDUCTION 概念図

MODEL ORDER REDUCTION 固有振動モードを利用

TTKT ξ

TT

m×1

m×n

m×n
n×n

F

g n ( A, B; u) = span {r0 , r1 , r2 ,L , rn −1}
n×1n×mK mor ξ = Fmor
K mor ξm×m m×1

 2次のKrylov部分空間を用いた直交化法（SOAR）

Fmorm×1

 以下の式でrrを求めていく。ただしC
Cは減衰マトリクス{r0 } = {K −1f K −1f }{r1} = [ A]{r0 }{ri } = [ A]{ri −1} + [B]{ri −2 }fori≥ 2

[ A ] = −[M ]−1[C][B] = −[M ]−1[K ]MODEL ORDER REDUCTION STATICMODEL ORDER REDUCTION STATIC 非線形FEMでの静的解析は

 MORは連立方程式を解く部分のオーダーを落とすことK (u)∆u = ∆F − Q(u)
 これにu = Tξ を代入し、左からT転置を掛けると

TT K (ξ )T∆ξ = TT (∆F − Q(ξ )) 整理すると次式

K mor ∆ξ = ∆Fmor − Q mor ただしK mor = TT K (ξ )T, ∆Fmor = TT ∆F, Q mor = TT Q(ξ )により計算時間を高速化する
 非線形有限要素解析では、剛性マトリクスの計算にかなりの時間を要する
 剛性マトリクスの計算を高速化する方法を検討

 ここで

K mor とQ mor を予め前処理としてなにかしらの近似関数で用意しておけば、K (ξ ) やQ(ξ ) を計算する
時間が省ける 応答曲面法の利用

MODEL ORDER REDUCTION STATICMODEL ORDER REDUCTION DYNAMIC

 強制変位解析

 動的解析

A∆ξ = ∆u
Π=&&t + Qt = F tMu

1 T∆ξ K mor ∆ξ − ∆ξT (∆Fmor − Q mor ) + λ T ( A∆ξ − ∆u)2δ Π = δ∆ξ (K mor ∆ξ − ∆Fmor + Q mor + A λ ) + δ λ ( A∆ξ − u )T K mor AA   ∆ξ   ∆Fmor − Q mor   = ∆u0  λ  T 強制変位解析

1L = ξ& T M mor ξ& + ξT (Fmor − Q mor ) + λ T ( Aξ − u )2M mor
 AAT   ξ t +∆t   2ξ t − ξ t −∆t + ∆t 2 (Fmor − Q mor )t =
2ut +∆t0   −∆t λ t  

中央差分

応答曲面法

応答曲面法

 計算時間短縮のため、

 いまk回計算を行なったとすると、線形回帰モデルは次予め応答曲面として用意できる部分を計算

K mor , Q mor 応答曲面とは、複数の変数とそこから予測される応答の関係式を近似したもの 出力がベクトルやマトリクスの場合でも、出力がスカラーの場合と同様の議論で応答曲面を作ることができる式で表される

 Q1   Ξ1  ε1  Q  Ξ  
 2  =  2  Β + ε 2  , Q = ΞΒ + ε M  M M    Q k  Ξk 
ε k 
 誤差の二乗和が最小になるようにBを決めると、

β 0 −1Β =  β1  = ( ΞT Ξ ) ΞT Q
β 2 

結果（精度）の評価

シミュレーションの検証 解析モデル

 FEMで一般的に用いられる誤差

 PCのCPUはIntel Core i7-2600 3.40GHz,メモリは8.0GB

変位

 図のようなモデルを用いる(軸の目盛りはmm){uFEM − u MOR } {uFEM − u MOR }E=T ヤング率0.05GPaT

ポアソン比0.3 密度930kg/m2{u FEM } {u FEM }

8080



エネルギー誤差707060{ε FEM − ε MOR } D {ε FEM − ε MOR }
E=60

T

50

{ε FEM }T D{ε FEM }50404030

 特徴を考慮した誤差

203010{u FEM − u MOR } [A]{u FEM − u MOR }E=[ A ] = diag (u FEM )T[ A ] = diag (u FEM )n{u FEM }T [ A]{u FEM }2002010200-20SOARによるMOR 強制変位解析

計算結果

 左：非線形FEM解析

 計算時間右：MOR808070706060505040403030202010100-30-20-1001020300-30100-10-200-30-20-10010

非線形FEM 170s、 MOR 0.1s
 通常の変位誤差
固有モード 荷重解析
固有モード 強制変位解析
SOAR
荷重解析
SOAR
強制変位解析

E0.5262
0.35190.56290.2949

 特徴を考慮した変位誤差-20-10010uFEM − uMOR
= 0.2949
uFEM

2030

固有モード 荷重解析
固有モード 強制変位解析
SOAR
荷重解析
SOAR
強制変位解析

E0.7205
0.15170.74860.13992030

シミュレーションの検証 動的解析

シミュレーションの検証 動的解析

 静的解析と同様の現象
 1秒間で線形的に強制変位を与える

8080707060606050505050404040403030303020202020101010100-300-30807060-20-100102030-20-1001020300-30-20-1001020300-30

150150100100505000-50

00.10.2-100100.3

0.40.5時間[s]0.6

nstep∆t∆tu iFEM− u iMORi =1∆tu iFEME=∑-2012345200ξξ 用いるモード数を5,10,20として検証

70250

12345200

 固有モードを用いたMORでの強制変位解析

80

MOR モード数5非線形FEM
2500.70.80.91-50

00.1

0.20.30.40.5時間[s]0.60.70.80.9nstep = 0.3560

20

動的解析 結果と考察

シミュレーションの検証 動的解析

 計算時間は、応答曲面を用いた場合モード数10では モード数の違いの検証を行なうために、支配的なモード0.22ｓほどであった
通常の変位誤差
モード数 E0.3560
50.3491
100.343920

特徴を考慮した変位誤差
モード数51020E0.1855
0.17870.1759が6番目以降に出てくる変形で検証
 椅子前面の1点に強制変位を垂直方向に与え、後方の1点に対して垂直方向を固定する強制変位を与えた
80807070606050405030204010

 少ないモードで入力を満たす必要があるため、支配的なモードが大きく出る傾向にある
 大まかな傾向は再現できている

1300402020402000-20-20-40-40100403020100-10-20-30-40

動的解析 リアルタイム性 検証

動的解析 リアルタイム性 考察

E[GPa
E[GPa]
GPa]

Δt[s]
Δt[s]

ステップ数

計算時間[s]
計算時間[s] 強制変位を与えたところで荷重をなくし、自由振動させ

0.05

0.001

1000

0.23400.10.001

1000

0.2496

0.20.001

1000

0.24960.3

0.00052000

0.48360.4

0.00052000

0.46800.8

0.00052000

0.4992るシミュレーション
 ヤング率1GPaでは振動数が大きすぎてまともなアニメーションにならなかった CGアニメーションとしての用途を考えると以下のアニメーションくらいが望まれると考えられる

1.00.0005

2000

0.4836

1.10.00025

4000

0.93600.001GPa

 1GPaのものまでリアルタイムで計算可能

ゴム0.001-0.1GPa テフロン0.5GPa ポリエチレン0.8GPa
ポリプロピレン1.5-2GPa ナイロン2-4GPa
出典：http://www.engineeringtoolbox.com

動的解析 リアルタイム性 考察

0.0001GPa よって十分堅さに関しては需要に応えうるリアルタイムシミュレーションが行なえていると言える

動的解析 リアルタイム性 考察

E[GPaE[GPa]
GPa]

Δt[s]
Δt[s]

ステップ数

計算時間[s]
計算時間[s]

0.001

0.01

1000.0312

 このヤング率ではさらに時間刻みを小さくとれる リアルタイムでの計算のボーダーラインの1/16程度で計算できる 連立方程式の解法はO(n3)であるので、単純に計算すれば理論上はモード数を25としてもリアルタイムで計算できる
 ただし今回の例(Δt=0.01)ではモード数23で発散してしまったため、22が限界のモード数となる

 応答曲面を作成するためにモード数を22とすると4^22回の計算が必要となり非常に計算コストがかかってしまう モード同士の相互作用を考慮し、直交表を用いるなどすることで計算回数を大幅に減らすことは可能である
 特に2次多項式で補間しているK morの応答曲面は複合計画と呼ばれる手法が非常に効果的である
 モード数を増やすことで、解の精度の向上と、安定性の向上が期待できる