📚

Material Point Method について

に公開

The Material Point Method Course

5 運動学の理論(KINEMATICS THEORY)

  • MPMで扱うパーティクルとは、物体の個々の粒子(原子や分子)ではなく、材料の連続体の中で局所的な性質を代表するサンプルである
  • 物質の各点で密度や速度、応力(内部で発生する力の状態)を連続的な関数として定義する。これにより、物体全体の変形や運動を微分方程式などを使って記述することができる

5-1 連続体運動(Continuum Motion)

  • 運動学(キネマティクス)とは、連続体(continuum)として扱われる物質における運動の研究を指す。ここでの主な焦点は、形状の変化、すなわち変形が局所的にも全体的にもどのように起こるかという点にある。

連続体力学では、変形は通常、「材料空間(あるいは未変形空間)\bm{X} 」と「世界空間(または変形後の空間)\bm{x}」、そしてこれらを結ぶ変形写像 \phi(X, t) を用いて表現される。ここで、 \bm{X} はシミュレーション開始時の各点の「初期位置」、\bm{x} は現在の「位置」と考えることができる。

\begin{equation} \bm{x} = \bm{x}(\bm{X}, t) = \phi(\bm{X}, t) \end{equation}

例えば、ある物体が一定の速度vで方向nに向かって移動している時は

\begin{equation} \bm{x} = \bm{X} + tv\bm{n} \end{equation}

また、物体が剛体(大きさを持ち、変形しない物体)であるとき、時刻 t における運動は

\begin{equation} \bm{x} = \bm{R}\bm{X} + \bm{b} \end{equation}

となる。剛体は並行(並進)運動と、回転運動を行う。\bm{R} は回転行列、\bm{b} は並進ベクトルである。

この変形写像 \phi は、連続体に基づく物理現象を定量化するために利用できる。例えば、ある材料点 \bm{X} の時刻 t における速度は

\begin{equation} \bm{V}(\bm{X}, t) = \frac{\partial \phi}{\partial t}(\bm{X}, t) \end{equation}

また、加速度は

\begin{equation} \bm{A}(\bm{X}, t) = \frac{\partial^2 \phi}{\partial t^2}(\bm{X}, t) = \frac{\partial \bm{V}}{\partial t}(\bm{X}, t) \end{equation}

となる。上記で定義された速度 \bm{V} や加速度 \bm{A} は、材料の位置 \bm{X} と時間 t の関数であり、これは「ラグランジュ視点」に基づいている。 ラグランジュ視点とは、物質の各点(材料点 \bm{X} )を固定して、その点に着目しながら、時間と共にその点がどのように動くかを観測することを意味する。これとは逆のものとして「オイラー視点」がある。オイラー視点では、空間内の固定された位置を基準に、そこを通過する材料展の性質(例えば速度)を測定するものである。

5-2 変形(deformation)

ここで、 \bm{X} を材料座標(未変形状態の座標)、 \bm{x} を世界座標(変形後の座標)とし、これらはそれぞれ領域 \Omega_{0}\Omega_{1} に属する。任意の材料点 \bm{x} \in \Omega_0 に対して、ある時刻 t\bm{x} = \phi(\bm{X}, t) という変形写像 \phi によって、 \bm{X}\Omega_t 上の対応する位置 \bm{x} に写される。

この変形写像 \phi のヤコビ行列は、以下で説明するいくつかの理由から非常に有用である。例えば、弾性体の物理は、このヤコビ行列を用いて自然に記述される。標準的な表記では、変形写像のヤコビ行列を \bm{F} と表し、

\begin{equation} \bm{F}(\bm{X}, t) = \frac{\partial \phi}{\partial \bm{X}}(\bm{X}, t) = \frac{\partial \bm{x}}{\partial \bm{X}}(\bm{X}, t) \end{equation}

と定義する。 \bm{F} はしばしば「変形勾配」と呼ばれ、離散的には 2 \times 2 (2次元の場合)または、 3 \times 3 (3次元の場合)の小さな行列となる。なお、3次元空間でも布や薄い殻状物体の場合、材料空間自体が実質2次元であるため、 \mathbf{F}3 \times 2 の行列となることもある。つまり、各材料点 \mathbf{X} に対して \mathbf{F}(\mathbf{X}, t) は、その点における局所的な典型の状態を記述する d \times d 行列( d = 2 または d = 3 )である。インデックス表示では、

\begin{equation} F_{i,j} = \frac{\partial \phi_i}{\partial X_i} = \frac{\partial x_i}{\partial X_i}, \quad i, j = 1, \dots , d \end{equation}

と表される。

例えば、式(2)の変形写像(物体が一定速度で平行移動する場合)の変形勾配は単位行列になる。また、式(3)の剛体運動(回転と並進のみ)の場合、変形勾配 \bm{F} は回転行列 \bm{R} となる。どちらの場合も、物体自体は変形しておらず、内部に応力が発生するべきではない。(ただし、芸術的な効果を狙う場合を除く)。

直感的には、変形勾配 \bm{F} は、材料が局所的にどの程度変形しているかを示している。例えば、シミュレーション開始時に材料内の2点 x^0_1x^0_2 があるとき、これらに対応する現在の位置が x_1x_2 であれば(図1)、

(x_2 - x_1)= \bm{F}(x^0_1 - x^0_2)

という関係が成立する。

また、 \bm{F} の行列式(determinant)、通称 J = det(\bm{F}) は非常に重要である。 J は、変形後の微小体積と初期の微小体積との比率を表しており、局所的な体積変化を特徴づける。例えば、剛体運動(回転・並進)の場合、 \bm{F} が回転行列となるため、 J = 1 となる。 J > 1 なら体積が増加し、 J < 1 なら体積が減少していることを意味する。

さらに、 J = 0 となると、体積がゼロになってしまうことを意味するが、実際の物理現象では起こり得ないものの、数値計算上はそのような \bm{F} が発生する可能性がある。3次元の場合、 J < 0 となると、物質が反転している状態(例えば、2次元の三角形で一つの頂点が対向する辺を突き抜け、面積が負になる)が起こる。こうした特異なケースに対しては、Invertible Elasticity [Irving et al., 2004; Stomakhin et al., 2012] のような手法が用いられることがある。

5-3 押し出しと引き戻し(Push Forward and Pull Back)

これまでの議論では、すべての量を材料座標 (\bm{X}, t) で表現してきた。これを「ラグランジュ視点」と呼ぶ。ここで用いる変形写像 \phi は全単射(bejective)であり、滑らかであると仮定しているため、初期領域 \Omega_0 と現在の領域 \Omega_t はお互いに位相同相(homeomorphic)、さらには微分同相(diffeomorphic)となる。すなわち、同じ時間において異なる材料点が同じ位置に来ることはなく、常に一位な対応が存在する。つまり、任意の \bm{x} \in \Omega_t に対して、必ず一つの \bm{X} \in \Omega_0 が存在し、 \phi (\bm{X}, t) = \bm{x} となるのである。これは逆写像 \phi^{-1}(·, t) : \Omega_t \rightarrow \Omega_0 が存在することを意味する。

この性質により、ある領域上に定義された関数は変数変換を通じて、別の領域上の関数として自然に考えることができる。これを「push forwad(押し出し)」または「pull back(引き戻し)」と呼ぶ。具体的には、 G : \Omega_0 \rightarrow \mathbb{R}g(\bm{x}, t) = G(\phi^{-1}(\bm{X}, t)) と定義される。同様に、 g の pull back は G(\bm{X}) = g(\phi(\bm{X}, t), t) と書くことができ、定義から G と一致する。

なお、 push forward された関数はしばしばオイラー視点の量( \bm{x} 座標で定義される)と呼ばれ、pull back された関数はラグランジュ視点の量( \bm{X} 座標で定義される)と呼ばれる。速度 \bm{V} (\bm{X}, t) や加速度 \bm{A}(\bm{X}, t) はラグランジュ視点で定義されてる。

\begin{gather} \bm{V}(\bm{X}, t) = \frac{\partial \phi}{\partial t}(\bm{X}, t) \\ \bm{A}(\bm{X}, t) = \frac{\partial^2 \phi}{\partial t^2}(\bm{X}, t) = \frac{\partial \bm{V}}{\partial t}(\bm{X}, t) \end{gather}

同様にオイラー視点の速度 \bm{v}(\bm{x}, t) や加速度 \bm{a}(\bm{x}, t) を次のように定義することもできる。

\begin{gather} \mathbf{v}(\bm{x}, t) = \bm{V}(\phi^{-1}(\bm{x}, t), t) \\ \bm{a}(\bm{x}, t) = \bm{A}(\phi^{-1}(\bm{x}, t), t) \end{gather}

また、 pull back の関係から、

\begin{gather} \bm{V}(\bm{X}, t) = \bm{v}(\phi(\bm{X}, t), t) \\ \bm{A}(\bm{X}, t) = \bm{a}(\phi(\bm{X}, t), t) \end{gather}

となる。さらに、連鎖律を用いると、ラグランジュ視点での加速度 \bm{A}(\bm{X}, t) は、

\begin{equation} \begin{split} \bm{A}(\bm{X}, t) &= \frac{\partial }{\partial t}\bm{V}(\bm{X}, t) = \frac{\partial \bm{v}}{\partial t}(\phi(\bm{X}, t), t) + \frac{\partial \bm{v}}{\partial \bm{x}}(\phi(\bm{X}, t), t)\frac{\partial \phi}{\partial t}(\bm{X}, t) \end{split} \end{equation}

と表され、インデックス記法では

\begin{equation} A_i(\bm{X}, t) = \frac{\partial}{\partial t}V_i(\bm{X}, t) = \frac{\partial v_i}{\partial t}(\phi(\bm{X}, t), t) + \Sigma_j \frac{\partial v_i}{\partial x_j}(\phi(\bm{X}, t), t)\frac{\partial \phi_j}{\partial t}(\bm{X}, t) \end{equation}

となる。

また、式(8)と式(10)を組み合わせると

\begin{equation} v_j(\bm{x}, t) = \frac{\partial \phi_j}{\partial t}(\phi^{-1}(\bm{x}, t), t) \end{equation}

となり、さらに式(11)と式(15)を組み合わせると、オイラー視点の加速度は

\begin{equation} a_i(\bm{x}, t) = A_i(\phi^{-1}(\bm{x}, t), t) = \frac{\partial v_i}{\partial t}(\bm{x}, t) + \frac{\partial v_i}{\partial x_j}(\bm{x}, t)v_j(\bm{x}, t) \end{equation}

となる。これらから、直感に反する以下の結果を得る

\begin{equation} a_i(\bm{x}, t) \neq \frac{\partial v_i}{\partial t}(\bm{x}, t) \end{equation}

5-4 物質微分(Material Derivative)

オイラー視点で定義された加速度 \bm{a}(\bm{x}, t) と速度 \bm{v}(\bm{x}, t) の関係は、単純な時間の偏微分だけでは表現できない。一方で、この関係は非常に一般的であり、「物質微分(material derivative)」と呼ばれる。物質微分の記法として、

\begin{equation} \frac{D}{Dt}v_i(\bm{x}, t) = \frac{\partial v_i}{\partial t}(\bm{x}, t) + \frac{\partial f}{\partial x_j}(\bm{x}, t)v_j(\bm{x}, t) \end{equation}

が導入され、これによって

\begin{equation} \bm{a} = \frac{D}{Dt}\bm{v} \end{equation}

となる。

一般のオイラー視点で定義された関数 f(·, t): \Omega^t \rightarrow \mathbb{R}
に対しても、同じ記法を用いて

\begin{equation} \frac{D}{Dt}f(\bm{x}, t) = \frac{\partial f}{\partial t}(\bm{x}, t) + \frac{\partial f}{\partial x_j}(\bm{x}, t)v_j(\bm{x}, t) \end{equation}

と定義する。ここで、 \frac{D}{Dt}f(\bm{x}, t) はラグランジュ視点で定義された関数 F(·,t): \Omega_0 \rightarrow \mathbb{R}
の時間微分 \frac{\partial F}{\partial t} のオイラー視点への押し出し(push forward)に相当する。つまり、 Ff の pull back であり、この \frac{\partial}{\partial t}F(\bm{X}, t)\frac{D}{Dt}f(\bm{x}, t)への押し出し規則は非常に有用で、常に念頭に置くべきある。

変形勾配 \bm{F} は通常ラグランジュ視点の量として考えられますが、オイラー視点への押し出し f としても扱うことができる。具体的には、 \bm{f}(·, t): \Omega_0 \rightarrow \mathbb{R}^{d \times d}\bm{F} の push forward としたとき、

\begin{equation} \frac{D}{Dt}\bm{f} = \frac{\partial \bm{v}}{\partial \bm{x}}\bm{f} \quad or \quad \frac{D}{Dt}f_{ij} = \Sigma_k \frac{\partial v_i}{\partial x_k}f_{kj} \end{equation}

となる。これは、流れ場の空間的な変化 \frac{\partial \bm{v}}{\partial \bm{x}} によって、 \bm{f} がどのように変化するかを示している。これが成り立つ理由は以下のように確認できる。

\begin{equation} \frac{\partial }{\partial t}F_{ij}(\bm{X}, t) = \frac{\partial }{\partial t}\frac{\partial \phi_i}{\partial X_j}(\bm{X}, t) = \frac{\partial V_i}{\partial X_j}(\bm{X}, t) = \frac{\partial v_i}{\partial x_k}(\phi(\bm{X}, t), t)\frac{\partial \phi_k}{\partial X_j}(\bm{X}, t) \end{equation}

これは、式(12)を微分することで得られる。

一部の文献([Bonet and Wood, 2008] や [Klar et al., 2016] など)では、上記の式(22)を \mathbf{f} の代わりに \bm{F} を用いて、

\begin{equation} \dot{\bm{F}} = (\nabla \bm{v})\bm{F} \quad or \quad \frac{D\bm{F}}{Dt} = (\nabla \bm{v})\bm{F} \end{equation}

と書く場合もある。また、 \dot{\bm{F}} = \frac{\partial}{\partial \bm{X}}(\frac{\partial \phi}{\partial t}) という形も現れる。このように書かれると、 \bm{F}\bm{f} の区別が曖昧になり、両空間での時間微分に \dot{\bm{F}} という記号が用いられる。文脈で \bm{F} = \bm{F}(\bm{X}, t) あるいは \bm{F} = \bm{F}(\bm{x}, t) と明確に示されていれば問題はないが、混乱を避けるため、通常はラグランジュ視点の量を \bm{F} 、オイラー視点の量を \bm{f} として区別することを好む。

最後に、式(23)は、各 MPM パーティクルにおける離散化された変形勾配 \bm{F} の更新式(Section 9-4 参照)を導出する上で非常に重要な役割を果たす。

5-5 体積と面積の変化(Volume and Area Change)

材料空間(未変形状態)に非常に小さな体積 dV が存在すると仮定する。ここで dV は標準基底ベクトル \bm{e}1, \bm{e}2, \bm{e}3 に沿って、微小な長さ dL_1, dL_2, dL_3 によって定義される体積要素は、 dV = dL_1\bm{e_1} \cdot (dL_2\bm{e_2} \times dL_3\bm{e_3}) と書ける(ただし、 dL_1, dL_2, dL_3 は非常に小さいスカラー量で、 dL_i = dL_ie_i と表される)。この場合、

\begin{equation} dV = dL_1 \cdot dL_2 \cdot dL_3 \end{equation}

この定義は、直交座標系での通常の微小体積の表現と一致する。

次に、材料領域にある各微小なベクトル d\bm{L}_i は、変形勾配 \bm{F} によって世界空間における対応する変形後の微小ベクトル d\bm{l}_i に写像される。

\begin{equation} d\bm{l}_i = F d\bm{L}_i \end{equation}

これらの変形後のベクトルから構成される体積は、 dl_1 dl_2 dl_3 = J dL_1 dL_2 dL_3 または dv = J dV となり、ここで J = det(\bm{F}) である。

この性質を利用すると、任意の関数 G(\bm{X}) または g(\bm{x}, t) に対し、変数変換を行うときに、押し出し(push forward)/ 引き戻し(pull back)の概念が良く使われる。すなわち、任意の領域 B^t \subseteq \Omega^t に対して、

\begin{equation} \int_{B_t} g(\bm{x}) d\bm{x} = \int_{B_0} G(\bm{X}) J(\bm{X}, t) d\bm{X} \end{equation}

と書ける。ここで、 B^0 は変形写像 \phi(\cdot, t) により B^t の前像(pre-image)であり、 Gg の引き戻し(pull back)である。 J(\bm{X}, t) は変形勾配の行列式である。

同様の議論は面積についても成り立つ。例えば、材料空間 \Omega^0 の任意の微小面積 d\bm{S} があり、その単位法線を \bm{N} とする。変形後の世界空間における対応する面積を d\bm{s} とし、その単位法線を \bm{n} とする。

\begin{equation} d\bm{S} = (dS) \bm{N} \end{equation}
\begin{equation} d\bm{s} = (ds) \bm{n} \end{equation}

これは、面積要素が大きさと法線方向の積で表現できることを示している。

さらに、面積 d\bm{S} ともう一つの非常に小さなベクトル d\bm{L} (その変形後の対応は d\bm{l} )を組み合わせると、体積は

\begin{equation} dV = d\bm{S} \cdot d\bm{L} \end{equation}
\begin{equation} dv = d\bm{s} \cdot d\bm{l} \end{equation}

ここで、内積を取ることにより、体積が得られる性質を用いています。

前述の結果 dv = J dV と組み合わせると、

\begin{equation} J d\bm{S} \cdot d\bm{L} = d\bm{s} \cdot (\bm{F} d\bm{L}) \end{equation}

ここで d\bm{l} = \bm{F} d\bm{L} を用いている。式(32)が任意の d\bm{L} に対して成立するため、次の関係が得られる。

\begin{equation} d\bm{s} = \bm{F}^{−T}Jd\bm{S} \end{equation}

または

\begin{equation} \bm{n}ds = \bm{F}^{-T}J\bm{N}dS \end{equation}

となります。この関係を用いることで、表面積分を以下のように記述できる。

\begin{equation} \int_{\partial B^t} \bm{h}(\bm{x}, t) \cdot \bm{n}(\bm{x})ds(\bm{x}) = \int_{\partial B^0} \bm{H}(\bm{X}) \cdot \bm{F}^{-T} (\bm{X}, t) \bm{N}(\bm{X})J(\bm{X}, t)dS(\bm{X}) \end{equation}

ここで、 \bm{H}:\Omega^0 \rightarrow \mathbb{R}^d\bm{h}:\Omega^t \rightarrow \mathbb{R}^d のpull back、 \bm{n}(\bm{x})\bm{x} における \partial B^t の単位外法線、 \bm{N}(\bm{X})\bm{X} における \partial B^0 の単位外法線です。これらの関係は、運動方程式を導出する際に非常に役立ちます。

6 超弾性(HYPERELASTICITY)

応力の物理的な意味を導入するのは、7章の運動方程式(運動量保存則)を導出する際の方が自然である。ここでは、応力はひずみ(この場合は変形勾配)と、何らかの「構成関係」を通じて関連しているという事実を紹介するにとどめる。応力は、領域全体に存在する場である。例えば、コーシー応力は \sigma(\cdot, t) : \Omega^t \rightarrow \mathbb{R}^{d \times d} と表される。離散的には、応力は評価された各点における小さなテンソル(行列)である。変形下での材料の応答を支配するためには、状態(変形勾配 \bm{F} など)と応力を関連付ける構成モデルが必要である。完全に超弾性的な材料の場合、構成関係はポテンシャルエネルギーを通じて定義され、初期状態からの非剛体変形とともに増加する。このセクションでは、弾性材料と、正確ではないものの実際には使いやすい塑性モデルに焦点を当てている。以下に説明するモデルは、MPMの論文や映画において、弾性体、砂、雪、溶岩、その他多くの材料のシミュレーションに成功裏に使用されている。

6-1 第一ピオラ-キルヒホッフ応力(First Piola-Kirchoff Stress)

伝統的な固体(弾性体など)を扱う際、ひずみと応力の関係を表すには、変形勾配(deformation gradient, \mathbf{F} )と第一ピオラ-キルヒホッフ応力(First Piola-Kirchoff stress, \mathbf{P} )を用いることが好まれる。なぜなら、これらの量は物体の変形前の初期状態(材料空間、material space)で自然に表現できるからである。超弾性体(Hyperelastic materials)とは、その第一ピオラ-キルヒホッフ応力 \mathbf{P} が、ひずみエネルギー密度関数(strain energy density function) \Psi (\bm{F}) から次のように導出できるような弾性体のことである。

\begin{equation} \bm{P} = \frac{\partial \Psi}{\partial \bm{F}} \end{equation}

添字表記を用いると、これは次のようになる。

\begin{equation} P_{ij} = \frac{\partial \psi}{\partial F_{ij}} \end{equation}

\psi(\bm{F}) は、非剛体的な(つまり形状や体積が変化する)変形勾配 \bm{F} に対してペナルティを与える(エネルギーが増加するように設計された)弾性エネルギー密度関数(スカラー関数)である。 \bm{P} は離散的には \bm{F} と同じ次元の小さな行列(テンソル)である。

工学分野でより一般的に用いられるコーシー応力(Cauchy stress, \mathbf{\sigma})と \bm{P} を関連付けるのは容易であり、次の式で表される。

\begin{equation} \bm{\sigma} = \frac{1}{J}\bm{P}\bm{F}^T = \frac{1}{det(\bm{F})} \frac{\partial \psi}{\partial \bm{F}}\bm{F}^T \end{equation}

材料の挙動は、変形写像 \phi と応力 \sigma または \bm{P} の相互作用によって定義される。超弾性体の場合、応力は変形勾配によって表現される形状の変化の関数である。材料の運動が剛体的(rigid)であるのは、次の場合である。

\begin{equation} \boldsymbol{\phi}(\bm{X}, t) = \bm{R}(t) \bm{X} + \bm{t}(t) \end{equation}

ここで、 $\bm{R}^T\bm{R} = \bm{I
}$ 、 det(\bm{R}) = 1 であり、 \bm{t} : [0, \infin) \rightarrow \mathbb{R}^d である。すなわち、 \bm{R} は回転、 \bm{t} は並進である。この場合、 \bm{F} = \bm{R} となることに注意が必要である。超弾性体は、 \bm{F} が直交行列(回転行列など)から逸脱することにペナルティを与えるエネルギーから生じる応力を介して、変形にペナルティを与えるのだ。これは次のようにして書くことができる。

\begin{equation} \bm{P} = \frac{\partial \Psi}{\partial \bm{F}}, \quad \Psi(\bm{F}) = \tilde{\Psi}(\bm{F}^T\bm{F}) \end{equation}

言い換えれば、 \bm{F} が直交行列であれば、エネルギーは変化せず(かつ最小値を持ち)、 \bm{F}^T\bm{F} はしばしば \bm{C}(右コーシー-グリーンテンソル, right Cauchy-Green tensor)で示される。もし材料が等方的(isotropic、つまり変形に対する応答が材料の方向に依存しない)であれば、エネルギーを \bm{C} の不変量(invariants)の関数として書くことで、さらに単純化できる。

\begin{equation} \Psi = \tilde{\Psi}(I_1, I_2, I_3) \end{equation}

ここで、 I_i\bm{F}^T\bm{F} の特性多項式の係数(しばしば等方性不変量と呼ばれる)であり、次のようになる。

\begin{gather} I_1 = tr(\bm{C}) \\ I_2 = tr(\bm{CC}) \\ I_3 = det(\bm{C}) = J^2 \end{gather}

コンピュータグラフィックス(CG)の分野では、これをさらに次のように書くのが便利である。

\begin{equation} \Psi(\bm{F}) = \hat{\Psi}(\boldsymbol{\Sigma}(\bm{F})) \end{equation}

ここで、 \bm{F} = \bm{U}\boldsymbol{\Sigma}\bm{Z}^T は、CGや力学分野で好まれる「Polar SVD(極特異値分解)」である([Irving et al., 2004; [source: 190] McAdams et al., 2011; Gast et al., 2016] 参照。なお、[McAdams et al., 2011] と [Gast et al., 2016] の両方で、これを高速に計算するためのオープンソースコードが公開されている)。これが「Polar SVD」と呼ばれるのは歴史的な理由からである。主に、 \bm{U}\bm{V} が回転行列であり、極分解(Polar decomposition) \bm{F} = \bm{R} \bm{S} が、 \bm{R} = \bm{U}\bm{V}^T および \bm{S} = \bm{V} \bm{Σ} \bm{V}^T を介して再構成できるからである。ここで、 \mathbf{R}\bm{F} に最も近い回転行列であり([Gast et al., 2016])、 \bm{S} は対称行列である。

等方性不変量 I_1I_2I_3 は特異値で書き表すことができるため、この方法は常に可能である。変形に対する構成応答(応力と変形の関係)をこのように構築することは、 \boldsymbol{\Sigma}(\bm{F}) を使うことで直感的に行うことができる。しかし、応力(および応力の微分)を得るためには、 \bm{F} の関数として特異値を微分する必要があり、これには注意深い導出が必要である。

6-2 ネオ・フック(Neo-Hookean)

ネオ・フック(Neo-Hookean)モデルは、弾性材料の大きな変形を予測するための、最も一般的な非線形超弾性モデルの一つである。このモデルのエネルギー密度関数は次のように与えられる。

\begin{equation} \Psi(\bm{F}) = \frac{\mu}{2}(tr(\bm{F}^T\bm{F}) - d) - \mu log(J) + \frac{\lambda}{2}log^2(J) \end{equation}

ここで、 d は問題の次元(2次元または3次元)を表し、 \mu\lambda はラメの定数と呼ばれ、ヤング率 E とポアソン比 \nu と次の関係がある。

\begin{equation} \mu = \frac{E}{2(1 + \nu)}, \quad \lambda = \frac{E \nu}{(1 + \nu)(1 - 2 \nu)} \end{equation}

\bm{F} が回転行列の場合、 \Psi(\bm{F}) = 0 となることは容易にわかる。 \bm{F} が反転していない(つまり J > 0)場合、 \Psi(\bm{F}) \geq 0 となる。このエネルギー密度関数は、超弾性固体を記述するのに適している。力の計算には、式(36)を用いて、 \bm{F} の関数として \bm{P} を得る必要があり、ネオ・フックモデルの場合、その結果は次のようになる。

\begin{equation} \bm{P} = \mu(\bm{F} - \bm{F}^{-T}) + \lambda log(J)\bm{F}^{-T} \end{equation}

陰解法による時間積分では、さらに \frac{\partial \bm{P}}{\partial \bm{F}} (応力の変形勾配による微分、すなわち応力テンソル)が必要になることに注意する。セクション6.4で、そのための実用的な方法を提供する。

6-3 固定コローテショナル構成モデル (Fixed Corotated Constitutive Model)

特異値分解(SVD)から定義される、もう一つの単純で広く使われているモデルがいわゆる固定コローテショナル(fixed corotated)モデルである。これが「固定(fixed)」と呼ばれる理由は、コンピュータグラフィックスの文献で一般的な、コローテショナル線形弾性(corotated linear elasticity)と呼ばれる広く使われているモデルに対する小さな修正版だからである。極特異値分解(polar SVD) \bm{F} = \bm{U} \bm{\Sigma} \bm{V}^{T} を仮定すると、固定コローテショナルモデルのエネルギーは次のようになる。

\begin{equation} \Psi(\bm{F}) = \hat{\Psi}(\bm{\Sigma}(\bm{F})) = \mu \sum_{i = 1}^{d}(\sigma_i - 1)^2 + \frac{\lambda}{2}(J - 1)^2 \end{equation}

ここで、もちろん J = \prod_{i=1}^{d} \sigma_i である。式中の \mu の項を展開すると、次のようになる。

\begin{equation} \mu \sum_{i=1}^{d}(\sigma_i - 1)^2 = \mu\left(\sum_{i=1}^{d} \sigma_i^2 - 2\sum_{i=1}^{d}\sigma_i + d\right) \end{equation}

そしてこれから、次のことを示せる。

\begin{equation} \frac{\partial}{\partial \bm{F}} \sum_{i=1}^{d}\sigma_i^2 = 2 \bm{F} \quad and \quad \frac{\partial}{\partial \bm{F}} \sum_{i=1}^{d}\sigma_i = \bm{R} \end{equation}

ここで、 \bm{F} = \bm{RS}\bm{F} の極分解(Polar decomposition)である( \bm{R} は回転行列、 \bm{S} は対称行列)。これはもちろん、 \bm{F} = \bm{U} \bm{V}^{T} \bm{V} \bm{\Sigma} \bm{V}^{T} のように \bm{F} の SVD から定義できる。すなわち、 \bm{R} = \bm{U} \bm{V}^T および \bm{S} = \bm{V} \bm{\Sigma} \bm{V}^T である。これらをすべて組み合わせると、次のようになる。

\begin{equation} \bm{P}(\bm{F}) = \frac{\partial \psi}{\partial \bm{F}}(\mathbf{F}) = 2\mu(\bm{F} - \bm{R}) + \lambda(J - 1)J\bm{F}^{-T} \end{equation}

二次微分(応力の微分)はもう少し注意が必要だが、比較的容易に計算できる。

まず、微分(differential、 \delta \bm{P} を結果として得る)を計算する。これは、陰解法ソルバーの行列フリー実装(matrix free implementation)に役立つ。

\begin{align} \frac{\partial^2 \Psi}{\partial \bm{F} \partial \bm{F}} : \delta \bm{F} &= \delta \left(\frac{\partial \Psi}{\partial \bm{F}} \right) \\ &= \delta(2 \mu (\bm{F} - \bm{R}) + \lambda(J - 1)J\bm{F}^{-T}) \\ &= 2\mu\delta\bm{F} - 2\mu\delta\bm{R} + \lambda J\bm{F}^{-T}\delta J + \lambda(J - 1)\delta(J\bm{F}^{-T}) \\ &= 2\mu\delta\bm{F} - 2\mu\delta\bm{R} + \lambda J\bm{F}^{-T}(J\bm{F}^{-T} : \delta\bm{F}) + \lambda(J - 1)\delta(J\bm{F}^{-T}) \end{align}

J\bm{F}^{-T} はその要素が \bm{F} の要素の多項式である行列なので、 \delta(J\bm{F}^{-T}) = \frac{\partial}{\partial \bm{F}}(J\bm{F}^{-T}) : \delta\bm{F} は直接容易に計算できる。残る課題は \delta \bm{R} の計算である。

\begin{align} \delta\bm{F} &= \delta\bm{R}\bm{S} + \bm{R}\delta\bm{S} \\ \bm{R}^T\delta\bm{F} &= (\bm{R}^T\delta\bm{R})\bm{S} + \delta\bm{S} \\ \bm{R}^T\delta\bm{F} - \delta\bm{F}^T\bm{R} &= (\bm{R}^T\delta\bm{R})\bm{S} + \bm{S}(\bm{R}^T\delta\bm{R}) \end{align}

ここでは、 \delta\bm{S} の対称性(Symmetry)と \delta\bm{R}^T\delta\bm{R} の歪対称性(Skew symmetry)を利用した。 \bm{R}^T\delta\bm{R} には3つの独立な成分があり、これらは直接解くことができる。方程式はこれらの成分に関して線形なので、 \bm{R}^T\delta\bm{R}3 \times 3 の連立方程式を解くことで計算できる。最後に、 \delta\bm{R} = \bm{R}(\bm{R}^T\delta\bm{R}) となる。

6-4 等方性弾性に対する実用的な微分戦略(A Practical Differentiation Strategy for Isotropic Elasticity)

ここでは、任意の一般的な等方性弾性材料に対して、第一ピオラ-キルヒホッフ応力 \bm{P} とその微分 \frac{\partial \bm{P}}{\partial \bm{F}} (テンソル形式または微分形式 \delta\bm{P} のいずれか)を計算するための実用的な方法([Stomakhin et al., 2012] に由来)を提供する。この方法は、記号計算ソフトウェアパッケージを利用する。ここでは Mathematica による実装について議論する。Maple や他のバージョンでも、同じロジックに従えば容易に作成できる。コンピュータグラフィックス応用で頻繁に発生する微分計算の実装に関するより詳細な議論については、[Schroeder, 2016] を参照すること。

この戦略は、 \frac{\partial \bm{P}}{\partial \bm{F}} と同様に見える他の対角空間での微分計算にも使用できることは注目に値する。例えば、いくつかのモデルでは、第一ピオラ-キルヒホッフ応力 \bm{P} の代わりにキルヒホッフ応力 \tau が使用される。

\begin{equation} \tau = \bm{U}\hat{\tau}\bm{U}^T \end{equation}

ここで、 \hat{\tau} は対角応力指標であり、各要素は \Sigma の関数である。 \frac{\partial \tau}{\partial \bm{F}} を計算するために、以下で説明する \frac{\partial \bm{P}}{\partial \bm{F}} の計算とほぼ全く同じ方法を使用できる。

6-4-1 \mathbf{P} の計算(Computing \mathbf{P}

まず \bm{P} から始める。等方性材料の場合、応力は次のように計算できる。

\begin{align} \bm{F} &= \bm{U}\bm{\Sigma}\bm{V}^T \\ \Psi(\bm{F}) &= \hat{\Psi}(\bm{\Sigma}) \\ \bm{P} &= \bm{U}\hat{\bm{P}}\bm{V}^T \end{align}

ここで、 \hat{\bm{P}} は対角行列であり、その要素は \frac{\partial \hat{\Psi}}{\partial \sigma_i} である。

以下は、対角空間での \bm{P} の証明(式63)である。まず、任意の回転 \bm{R} に対して \bm{P}(\bm{RF}) = \bm{RP}(\bm{F}) を示したい。任意のモデル(等方性である必要さえない)を考える。変形後の回転はエネルギーを変えるべきではない。これは \Psi(\bm{F}) = \Psi(\bm{RF}) を意味する。この方程式の微分を取ると、次のようになる。

\begin{align*} \delta\Psi = \frac{\partial \Psi}{\partial \bm{F}}(\bm{F}) : \delta\bm{F} &= \frac{\partial \Psi}{\partial \bm{F}}(\bm{RF}) : \delta(\bm{RF}) \\ (\bm{P}(\bm{F})) : (\delta\bm{F}) &= (\bm{P}(\bm{RF})) : \delta(\bm{RF}) \\ (\bm{P}(\bm{F})) : (\delta\bm{F}) &= (\bm{P}(\bm{RF}))_{ij} R_{ik} \delta F_{kj} \\ (\bm{P}(\bm{F})) : (\delta\bm{F}) &= (\bm{R}^T \bm{P}(\bm{RF})) : \delta(\bm{F}) \\ \bm{P}(\bm{F}) &= \bm{R}^T\bm{P}(\bm{RF}) \\ \bm{RP}(\bm{F}) &= \bm{P}(\bm{RF}) \end{align*}

同様に、等方性材料を仮定すれば、 \Psi(\bm{FR}) = \Psi(\bm{F}) を用いて、任意の回転 \bm{R} に対して \bm{P}(\bm{FR}) = \bm{P}(\bm{F})\bm{R} を証明できる。 \mathbf{P} に対するこれら二つの性質を用いると、次のようになる。

\begin{equation} \bm{P}(\bm{F}) = \bm{P}(\bm{U} \Sigma \bm{V}^T) = \bm{UP}(\Sigma)\bm{V}^T = \bm{U}\hat{\bm{P}}\bm{V}^T \end{equation}

6-4-2 \frac{\partial \bm{P}}{\partial \bm{F}} または \delta\bm{P} の計算 (Computing \frac{\partial \bm{P}}{\partial \bm{F}} or \delta\mathbf{P})

アイデアは \bm{P} の場合と同様であり、任意の二つの回転 \bm{R}\bm{Q} を取り、 \bm{P} の結果を用いると、次のようになる。

\begin{equation} \bm{P}(\bm{F}) = \bm{P}(\bm{R}\bm{R}^T\bm{F}\bm{Q}\bm{Q}^T) = \bm{R}\bm{P}(\bm{R}^T\bm{F}\bm{Q})\bm{Q}^T \end{equation}

\bm{K} = \bm{R}^T\bm{FQ} をすれば、以下を得る。

\begin{equation} \bm{P}(\bm{F}) = \bm{R}\bm{P}(\bm{K})\bm{Q}^T \end{equation}

\bm{R}\bm{Q} は定数であることに注意し、微分をとると、

\begin{align} \delta\bm{P} &= \bm{R} \left[\frac{\partial \bm{P}}{\partial \bm{F}}(\bm{K}) : \delta(\bm{K}) \right]\bm{Q}^T \\ &= \bm{R} \left[\frac{\partial \bm{P}}{\partial \bm{F}}(\bm{K}) : (\bm{R}^T\delta\bm{FQ}) \right]\bm{Q}^T \end{align}

\bm{R}\bm{Q} は自由に選べるため、 \bm{R} = \bm{U} かつ \bm{Q} = \bm{V} と選ぶ。すると \bm{K} = \Sigma となる。式は次のようになる。

\begin{equation} \delta\bm{P} = \bm{U}\left[\frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) : (\bm{U}^T\delta\bm{FV}) \right]\bm{V}^T \end{equation}

これは \delta\bm{P} を与えます。テンソル \frac{\partial \bm{P}}{\partial \bm{F}} については、添字表記を採用する。

\begin{align} (\delta\bm{P}_{ij}) &= U_{ik}\left(\frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) \right)_{klmn}U_{rm}\delta F_{rs}V_{sn}V_{jl} \\ (\delta\bm{P}_{ij}) &= \left(\frac{\partial \bm{P}}{\partial \bm{F}}(\bm{F}) \right)_{ijrs}\delta F_{rs} \end{align}

これら二つの式は任意の \delta\bm{F} に対して成り立たなければならないので、次のことが明らかになる。

\begin{equation} \left(\frac{\partial \bm{P}}{\partial \bm{F}}(\bm{F}) \right)_{ijrs} = \left(\frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) \right)_{klmn}U_{ik}U_{rm}V_{sn}V _{il} \end{equation}

したがって、残る問題は \frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) を計算することである。これを3次元で行う方法を示す。

まず、ロドリゲスの回転公式を使う必要がある。それは、任意の回転行列が単位ベクトル \bm{k} と回転角 \theta で書けることを示している。

\begin{equation} \bm{R} = \bm{I} + \sin(\theta)\bm{K} + (1 - \cos(\theta))\bm{K}^2 \end{equation}

ここで、 \bm{K}\bm{k} の歪対称外積行列である。これは、すべての回転行列がわずか3つの自由度 r1, r2, r3 しか持たないことを意味する。ここで、 \bm{k} = \frac{\bm{r}}{|\bm{r}|} かつ \theta = |\bm{r}| である。したがって、回転行列 \bm{U}\bm{V} をそれぞれ3つの数値でパラメータ化することができる。

これで、 s1,s2,s3,u1,u2,u3,v1,v2,v3  で \bm{F} を定義するための以下のコードがあります。ここで、 \bm{U}\bm{V} はロドリゲスの回転公式を用いて u_iv_i で定義され、 s_i\Sigma からの特異値である。

(Mathematicaコード例:3D)

id = IdentityMatrix[3];
var = {s1, s2, s3, u1, u2, u3, v1, v2, v3};
Sigma = DiagonalMatrix[{s1 ,s2 ,s3}];
cp[k1_, k2_, k3_] = {{0, -k3, k2}, {k3, 0, -k1}, {-k2, k1, 0}};
vV = {v1, v2, v3};
vU = {u1, u2, u3};
nv =  Sqrt[Dot[vV, vV]];
nu = Sqrt[Dot[vU, vU]];
UU = cp[u1, u2, u3] / nu;
VV = cp[v1 ,v2 ,v3] / nv;
U = id + Sin[nu] * UU + (1 - Cos[nu]) * UU.UU;
V = id + Sin[nv] * VV + (1 - Cos[nv]) * VV.VV;
F = U.Sigma.Transpose[V];

ここで、 cp は外積行列を生成するための関数である(式(73)の \bm{K} の計算に相当)。

これ以降、 3 \times 3 \times 3 \times 3 テンソル \frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) や他のそのようなテンソルを 9 \times 9 行列として記述していく。これは、各 3 \times 3 行列がサイズ9のベクトルになることを意味する。以前の \frac{\partial P_{ij}}{\partial F_{kl}}\frac{\partial P_{3(i-1)+j}}{\partial F_{3(k-1)+l}} になることは容易にわかる。さらに、ベクトル \bm{S} ={s1,s2,s3,u1,u2,u3,v1,v2,v3}\bm{F} のパラメータ化と呼ぶ。すると、連鎖律を適用することができ、

\begin{equation} \frac{\partial \bm{P}}{\partial \bm{F}}(\Sigma) = \frac{\partial \bm{P}}{\partial \bm{S}}(\Sigma)\frac{\partial \bm{S}}{\partial \bm{F}}(\Sigma) \end{equation}

となる。以下はこれを計算したコードの例である。ほぼ0回転に相当する極限 { u1, u2, u3, v1, v2, v3 } = + \epsilon をとることによって、 \bm{F} = \Sigma を達成していることに注意する。

(Mathematicaコード例:微分計算と連鎖律)

dFdS = D[Flatten[F], {var}];
dFdS0 = dFdS / .{u1->e, u2->e, u3->e, v1->e, v2->e, v3->e};
dFdS1 = Limit[dFdS0, e->0, Direction->-1];
dSdF0 = Inverse[dFdS1];
Phat = DiagonalMatrix[{t1[s1, s2, s3], t2[s1, s2, s3], t3[s1, s2, s3]}];
P = U.Phat.Transpose[V];
dPdS = D[Flatten[P], {var}];
dPdS0 = dPdS / .{u1->e, u2->e, u3->e, v1->e, v2->e, v3->e};
dPdS1 = Limit[dPdS0, e->0, Direction->-1];
dPdF = Simplify[dPdS1.dSdF0];

Mathematicaにおける”Direction → -1” は大きな値から小さな極限値へと近づく極限をとることを意味する。Mathematica による計算結果は、特異値と \hat{\bm{P}}
(対角応力) の項で与えられる。その後、その数式をコードに実装するために利用できる。[Stomakhin et al., 2012] は、 \frac{\partial \bm{P}}{\partial \bm{F}} (サイズ 9 \times 9 行列)が並べ替えられて、対角ブロック $\bm{A}^{3\times3}, \bm{B}{12}^{2\times2}, \bm{B}{13}^{2\times2},
\bm{B}_{23}^{2\times2}$ を持つブロック対角行列になる結果を与える。これは

\begin{equation} \bm{A} = \begin{pmatrix} \hat{\psi}_{,\sigma_1\sigma_1} & \hat{\psi}_{,\sigma_1\sigma_2} & \hat{\psi}_{,\sigma_1\sigma_3} \\ \hat{\psi}_{,\sigma_2\sigma_1} & \hat{\psi}_{,\sigma_2\sigma_2} & \hat{\psi}_{,\sigma_2\sigma_3} \\ \hat{\psi}_{,\sigma_3\sigma_1} & \hat{\psi}_{,\sigma_3\sigma_2} & \hat{\psi}_{,\sigma_3\sigma_3} \end{pmatrix} \end{equation}

\begin{equation} \bm{B}_{ij} = \frac{1}{\sigma_i^2 - \sigma_j^2}\begin{pmatrix}\sigma_i\hat{\psi}_{,\sigma_i}-\sigma_j\hat{\psi}_{,\sigma_j} & \sigma_j\hat{\psi}_{,\sigma_i} -\sigma_i\hat{\psi}_{,\sigma_j} \\ \sigma_j\hat{\psi}_{,\sigma_i} -\sigma_i\hat{\psi}_{,\sigma_j} & \sigma_i\hat{\psi}_{,\sigma_i} -\sigma_j\hat{\psi}_{,\sigma_j} \end{pmatrix} \end{equation}

になる。 \sigma_i^2 - \sigma_j^2 による除算は、二つの特異値がほぼ等しい場合、または二つの特異値の和がほぼゼロになる場合に問題となる。後者は、負の特異値を許容する慣例(可逆弾性 [Irving et al., 2004; Stomakhin et al., 2012] など)で可能である。

\bm{B}_{ij} を部分分数に展開すると、有用な分解が得られる。

\begin{equation} \bm{B}_{ij} = \frac{1}{2}\frac{\hat{\psi}_{,\sigma_i} - \hat{\psi}_{,\sigma_j}}{\sigma_i - \sigma_j} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} + \frac{1}{2}\frac{\hat{\psi}_{,\sigma_i} + \hat{\psi}_{,\sigma_j}}{\sigma_i + \sigma_j} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix} \end{equation}

\hat{\psi} が特異値の置換に対して不変であれば、 \sigma_i \rightarrow \sigma_j のとき \hat{\psi}_{,\sigma_i} \rightarrow \hat{\psi}_{,\sigma_j} となることに注意する。したがって、第一項は通常、注意深く実装すれば等方性モデルに対してロバストに計算することができる。もう一方の分数は、より深い意味合いを持つ。この項は、 \sigma_i + \sigma_j \rightarrow 0 のときに \hat{\psi}_{\sigma_i} + \hat{\psi}_{,\sigma_j} \rightarrow 0 であればロバストに計算ができる。この特性は好ましくはなく、なぜなら、構成モデルが多くの反転した構成から回復するのが困難になることを意味するからである。我々は通常、反転に対してロバストな挙動を示すモデルに興味があるので、この項は必然的にいくつかの状況下で非有界になる。これに対処するために、除算の前に分母の絶対値が 10^{-6} より小さくならないようにクランプし、微分を制限する。

2次元の場合、回転行列は単一の \theta (パラメータ) で単純にパラメータ化され、その再構成は次のようになる。

\bm{R} = \begin{pmatrix} \cos{\theta} & -\sin{\theta} \\ \sin{\theta} & \cos{\theta} \end{pmatrix}

2次元の場合のMathematicaコードは以下である。

id = IdentityMatrix[2];
var = {s1 ,s2 ,u1 ,v1};
S = DiagonalMatrix[{s1, s2}];
U = {{Cos[u1], -Sin[u1]}, {Sin[u1], Cos[u1]}};
V = {{Cos[v1], -Sin[v1]}, {Sin[v1], Cos[v1]}};
F = U.S.Transpose[V];
dFdS = D[Flatten[F], {var}];
dFdS0 = dFdS /. {u1->e, v1->e};
dFdS1 = Limit[dFdS0, e->0, Direction->-1];
dSdF0 = Inverse[dFdS1];
Phat = DiagonalMatrix[{t1[s1, s2], t2[s1 ,s2]}];
P = U.Phat.Transpose[V];
dPdS = D[Flatten[P], {var}];
dPdS0 = dPdS /. {u1->e ,v1->e};
dPdS1 = Limit[dPdS0, e->0, Direction->-1];
dPdF = Simplify[dPdS1.dSdF0];

6-5 雪の塑性(Snow Plasticity)

雪の構成挙動(力に対する応答)は、非常に広範で複雑な要因に依存する。対象とする条件に応じて、多くの異なるモデルが使用される。動的な挙動の多くは、比較的単純な弾塑性(elasto-plastic)の仮定で近似することができる。ここでは、硬化(hardening)効果と組み合わせた、単純な大変形塑性流動モデルを使用することにする。

塑性を表現するために、変形勾配 \bm{F} を弾性部分 \bm{F}_E と塑性部分 \bm{F}_P に次のように分解する。

\begin{equation} \bm{F} = \bm{F}_E\bm{F}_P \end{equation}

変形勾配は、材料がその運動によって局所的にどれだけ回転し変形したかの尺度である。変形勾配をこのように分解することで、この変形履歴を二つの部分に分ける。塑性部分 \bm{F}_P は、材料の履歴のうち「忘れ去られた」部分を表している。もし金属棒を曲げてコイル状のバネにした場合、棒は元々まっすぐだったことを忘れていると言え、そのコイル状のバネは、常にコイル状であったかのように振る舞う。この操作に関わるねじれや曲げは \bm{F}_P に保存される。もしバネがわずかに圧縮されると、バネはひずみ(変形)を感じる。これが弾性変形であり、 \bm{F}_E に保存される。バネはこの変形を覚えており、応答として材料はコイル状の形状に戻ろうとして応力を発生させる。このようにして、応力の計算には \bm{F}_E のみを使用すべきであることがわかる。金属棒の全履歴は、バネの形状に曲げられ (\bm{F}_P) 、その後圧縮される(\bm{F}_E)ことで構成される。

弾性応答は \bm{F}_E のみの関数である。直感的には、これは \bm{F}_P への局所的な移行における変形が永続的であることを示している。ある意味で、 \bm{F}_P は材料の新しい局所的な静止状態を表すようになる。この永続的な変形への移行は、通常、大きな変形に応じて起こる。この簡単な例は、アルミ缶をへこませることであろう。一度へこむと、すべての弾性応答は、へこんだ形状からの変位に対するものになる。ほとんどのモデルでは、応力に基づく条件に応じてこの分解を定義する。しかし、我々の場合、制約を \bm{F} (訳注: \bm{F}_E の誤りと思われる) 自体の特異値で定義されると考える方がより直感的であろう。これにより、塑性効果に対する視覚的な制御が向上する。具体的には、 \bm{F}_E の特異値 \sigma_{Ei} が、いくつかの小さな定数 \theta_c\theta_s に対して、区間 [1 - \theta_c, 1 + \theta_s] 内にあるように強制する。これは以下の手順で行われます。ある時刻 n

\begin{equation} \bm{F}^n = \bm{F}_E^n\bm{F}_P^n \end{equation}

が与えられ(ここで \bm{F}_E^n の特異値 \sigma_{Ei}^n は区間 [1 - \theta_c, 1 + \theta_s] 内にある制約を満たす)、新しい変形勾配 \bm{F}^{n+1} があるとする。まず、 \bm{F}^n から \bm{F}^{n+1}
への遷移で導入された新しい変形はすべて弾性的であったと仮定する。すなわち、まず新しい \bm{F}^{n+1} が次のように分解できると仮定する。

\begin{equation} \bm{F}^{n+1} = \hat{\bm{F}}_E^{n+1}\bm{F}_P^n \end{equation}

そうすることで、これは \hat{\bm{F}}_E^{n+1} を次のように定義する。

\begin{equation} \hat{\bm{F}}^{n+1} = \bm{F}_E^{n+1}(\bm{F}_P^n)^{-1} \end{equation}

実際には、完全な \bm{F} を保存する代わりに \bm{F}_E\bm{F}_P を保存する方が便利なことが多い。 \bm{F} の仮の更新は、 \bm{F}_E に直接適用して \hat{\bm{F}}_E^{n+1} を与えることができる。次のステップは、 \hat{\bm{F}}_E^{n+1} の特異値が区間 [1 - \theta_c, 1 + \theta_s] 内にあるという制約を強制することである。すなわち、\hat{\bm{F}}_E^{n+1} の特異値を \hat{\sigma}_{Ei}^{n+1}として、次のように定義する。

\begin{equation} \sigma_{E_i}^{n+1} = clamp(\hat{\sigma}_{E_i}^{n+1}, 1-\theta_c, 1+\theta_s), \quad i = 1, ...,d \end{equation}

ここで、 \hat{\bm{F}}_E^{n+1} の特異値分解が次であると仮定すると、

\begin{equation} \hat{\bm{F}}_E^{n+1} = \bm{U}_E^{n+1}\hat{\Sigma}_E{\bm{V}_E^{n+1}}^T \end{equation}

クランプされた特異値 \Sigma_E^{n+1} から \bm{F}_E^{n+1} を次のように定義できる。

\begin{equation} \bm{F}_E^{n+1} = \bm{U}_E^{n+1}\Sigma_E^{n+1}{\bm{V}_E^{n+1}}^T \end{equation}

もちろん、この \bm{F}_E^{n+1} の定義では、 \bm{F} の同じ分解を維持する必要があるため、次を満たすような新しい \bm{F}_P^{n+1} を決定する必要がある。

\begin{equation} \bm{F}^{n+1} = \bm{F}_E^{n+1}\bm{F}_P^{n+1} \end{equation}

しかしもちろん、 \bm{F}^{n+1} と新しい \bm{F}_E^{n+1} がわかっていれば、変形勾配の塑性成分は次のようになる。

\begin{equation} \bm{F}_P^{n+1} = \left(\bm{F}_E^{n+1}\right)^{-1}\bm{F}^{n+1} \end{equation}

雪は圧縮されるとより硬くなる傾向がある。の現象はしばしば硬化(hardening)と呼ばれ、我々はこの効果を加えるために、構成モデルに簡単な修正を加える。具体的には、ラメ係数 \mu\lambda が圧縮下で増加し、伸張下で減少するようにする。伸張下での材料強度の低下は、雪の分裂や破壊を促進する。これは広範な視覚現象にとって重要な特性です。この硬化効果を次のように定量化する。

\begin{equation} \mu(\bm{F}_P) = \mu_0e^{\xi(1 - J_P)}, \quad \lambda(\bm{F}_P)=\lambda_0e^{\xi(1-J_P) } \end{equation}

ここで \mu_0\lambda_0 は、元のヤング率とポアソン比から設定されたラメパラメータである。 \xi は硬化パラメータであり、通常 3 から 10 の間の値を使用する。また、 J_P = \det(\bm{P}) である。 \mu\lambda が大きくなりすぎないように、何らかの安全策を含めることも重要であり、これは、 \xi(1 - J_P) に数値的なクランプ上限を設けることで簡単に行うことができる。

また、導出からわかるように、 \bm{F}_P を保存する代わりに、 J_P を追跡するだけで十分である。ただし、時刻 n の量から \hat{\bm{F}}_E^{n+1} を取得する方法を知っている必要がある。MPM では、これは実際に可能である。

この硬化モデルは、雪が引き伸ばされると柔らかくなり(破壊を可能にするため)、圧縮されると雪玉を固めるように硬くなる という直感に基づいて設計されている。式(87)のルールは単なる経験的な式である。ルールを創造的にすることで、より多様で芸術的な材料の挙動を生み出すのに役立つ。

より厳密に導出された塑性モデルは、熱力学の第二法則に従う必要があり、また J_P=1(塑性変形は材料の体積を変えない)を強制する必要がある。この問題に関するより詳細な議論については、[Klar et al., 2016] の技術文書を参照すること。

7 支配方程式(GOVERNING EQUATIONS)

対象とする支配方程式は、質量保存則と運動量保存則である。以下にその結果を示し、その導出をセクション7.1と7.2で提供する。さらに、力のつり合いの弱形式(weak form)をセクション7.3で導出する。弱形式は、セクション9で方程式の最終的な時間的および空間的な離散化を導出するために不可欠である。導出を見る前に、セクション5(Kinematics Theory)を復習することを勧める。

\bm{X} 上で定義された速度を \bm{V}(\bm{X}, t) = \frac{\partial \Phi(\bm{X}, t)}{\partial t}=\frac{\partial \bm{x}(\bm{X}, t)}{\partial t} とすると、方程式のラグランジュ表示(Lagrangian view)は次のようになる [Gonzalez and Stuart, 2008]。

\begin{gather} R(\bm{X}, t)J(\bm{X}, t) = R(\bm{X}, 0) \quad \text{Conservation of mass} \\ R(\bm{X}, 0)\frac{\partial \bm{V}}{\partial t} = \nabla^T\cdot\bm{P}+R(\bm{X}, 0)\bm{g} \quad \text{Conservation of momentum} \end{gather}

ここで、 \bm{X} \in \Omega_0, t > 0 である。ここで R はラグランジュ質量密度であり、より一般的に使用されるオイラー質量密度 \rho とは R(\bm{X}, t) = \rho(\Phi(\bm{X}, t), t) の関係がある。質量保存則は \frac{\partial}{\partial t}(R(\bm{X}, t)J(\bm{X}, t)) = 0 とも書けることに注意すること。

オイラー表示(Eulerian view)では、支配方程式は次のようになる。

\begin{gather} \frac{D}{Dt}\rho(\bm{x}, t) + \rho(\bm{x}, t)\nabla^x\cdot\bm{v}(\bm{x}, t) = 0 \quad \text{Conservation of mass} \\ \rho(\bm{x}, t)\frac{D\bm{v}}{Dt} = \nabla^x\cdot\boldsymbol{\sigma} + \rho(\bm{x}, t)\bm{g} \quad \text{Conservation of momentum} \end{gather}

ここで、 \bm{v} = \bm{v}(\bm{x}, t) はオイラー速度(Eulerian velocity)、 \frac{D}{Dt}=\frac{\partial}{\partial t}+\bm{v}\cdot\nabla^x は物質微分(material derivative)である。

運動量保存則で示されるように、ラグランジュ表示では第一ピオラ-キルヒホッフ応力 \bm{P} を、オイラー表示ではコーシー応力 \bm{\sigma} を使用する方が便利である。

7-1 質量保存則(Conservation of Mass)

\rho(\bm{x}, t) をオイラー質量密度とし、 R(\bm{X}, t) をそのラグランジュ対応物(プルバック)とすると、次の関係がある。

\begin{gather} R(\bm{X}, t) = \rho(\Phi(\bm{X}, t), t) \\ \rho(\bm{x}, t) = R(\Phi^{-1}(\bm{x}, t), t) \end{gather}

\rho\Omega^t 上で自然に次のように定義されると考えることができる。

\begin{equation} \rho(\bm{x}, t) = \lim_{\epsilon \rightarrow +0}\frac{\text{mass}(B_\epsilon^t)}{\int_{B_\epsilon^t}\mathrm{d}\bm{x}} \end{equation}

ここで B_\epsilon^t は任意の \bm{x} \in \Omega^t を囲む半径 \epsilon の球である。 \text{mass}(B_\epsilon^t) は測定可能な量であるべきなので、これは間違いなく自然な定義である。質量保存則は次のように表現できる。

\begin{align} \text{mass}(B_\epsilon^t) &= \int_{B_\epsilon^t}\rho(\bm{x}, t)\mathrm{d}\bm{x} \nonumber\\ &= \int_{B_\epsilon^t}R(\bm{X}, t)J\mathrm{d}\bm{X} \nonumber\\ &= \int_{B_\epsilon^t}R(\bm{X}, 0)\mathrm{d}\bm{X} \nonumber\\ &= \text{mass}(B_\epsilon^0) \end{align}

すべての B_\epsilon^t \subset \Omega^t に対して(以前と同様、 B_\epsilon^0\phi(\cdot, t) の下での B_\epsilon^t の原像である)。二番目の等式は、セクション5.5の変数変換の結果を使用することから来ている。

これは単に、 B_\epsilon^t 内の質量(質量密度の積分として表現される)は時間とともに変化しないはずだと言っている。この集合は時刻 t における材料の一部に関連付けられており、流れの中で進化するにつれて、材料はより多くの、あるいはより少ない空間を占めるようになりますが、集合内の材料が多くなったり少なくなったりすることはないのである。 B_\epsilon^t は任意なので、次が真でなければならない。

\begin{equation} R(\bm{X}, t)J(\bm{X}, t) = R(\bm{X}, 0), \quad \forall{\bm{X}} \in \Omega^0, \quad t \geq 0 \end{equation}

J(\bm{X}, 0) = 1 であることに注意する。代わりに、質量保存則は次のように書くこともできる。

\begin{equation} \frac{\partial}{\partial t}(R(\bm{X}, t)J(\bm{X}, t)) = 0 \end{equation}

オイラー表示に切り替えるには、まず次のことに気づく。

\begin{equation} \frac{\partial}{\partial t}(RJ) = \frac{\partial R}{\partial t}J + R\frac{\partial J}{\partial t} \end{equation}

また、

\begin{equation} \frac{\partial J}{\partial t} = \frac{\partial J}{\partial F_{ij}}\frac{\partial F_{ij}}{\partial t} = JF_{ji}^{-1}\frac{\partial V_i}{\partial X_j} = JF_{ji}^{-1}\frac{\partial v_i}{\partial x_k}F_{kj} = J\delta_{ik}\frac{\partial v_i}{\partial x_k} = J\frac{\partial v_i}{\partial x_i} \end{equation}

ここで、行列式の微分規則を使用した。

\begin{equation} \frac{\partial J}{\partial \bm{F}} = J\bm{F}^{-T} \end{equation}

式(97)、(98)、(99) を組み合わせると、次のようになる。

\begin{equation} \frac{\partial R}{\partial t}J + Rj\frac{\partial v_i}{\partial x_i} = 0 \end{equation}

両辺をプッシュフォワード(push forward)すると、質量保存則のオイラー表示が得られる。

\begin{equation} \frac{D}{Dt}\rho(\bm{x}, t) + \rho(\bm{x}, t)\nabla \cdot\bm{v}(\bm{x}, t) = 0, \quad \forall\bm{x} \in \Omega^t, \quad t \geqq 0 \end{equation}

7-2 運動量保存則(Conservation of Momentum)

連続体に働く力は、体積力(body forces、例: 重力)または表面力(surface forces、応力ベース)のいずれかに分類される。応力ベースの力は、まずトラクション場(traction field)を介して定義され、その存在を仮定する。単位面積あたりの力(またはトラクション)場 \bm{t}(\cdot, \bm{n}, t) : \Omega^t \rightarrow \mathbb{R}^d は、次の関係によって定義される。

\begin{equation} \text{force}_s(B_{\epsilon}^t) = \int_{\partial B_\epsilon^t}\bm{t}(\bm{x}, \bm{n}(\bm{x}))\mathrm{d}s(\bm{x}) \end{equation}

ここで \text{force}_s(B_\epsilon^t) は、任意の領域 B_\epsilon^t に対して、その境界 \partial B_\epsilon^t の外側の材料から内側の材料へ及ぼされる正味の力である。すなわち、 \bm{t}(\bm{x}, \bm{n}, t) は、点 \bm{x} における材料の \bm{n}^+ 側(法線 \bm{n} の方向)の材料が \bm{n}^- 側の材料に及ぼす単位面積(d = 3)/ 単位長さ(d = 2)あたりの力である。これは、次を満たす応力場(コーシー応力, Cauchy stress) \sigma(\cdot, t) : \Omega^t \rightarrow \mathbb{R}^{d \times d} の存在を意味することが示せる。

\begin{equation} \bm{t}(\bm{x}, \bm{n}, t) = \sigma(\bm{x}, t) \bm{n} \end{equation}

\bm{v} をオイラー速度(ラグランジュ対応物は \bm{V})とする。すると、運動量保存則は次のように表現される。

\begin{equation} \begin{aligned} \frac{d}{dt}\int_{B_\epsilon^t}\rho(\bm{x}, t)\bm{v}(\bm{x}, t)d\bm{x} &= \frac{d}{dt}\int_{B_\epsilon^t}R(\bm{X}, t)\bm{V}(\bm{X}, t)J\mathrm{d}\bm{X} \\ &= \int_{B_\epsilon^0}R(\bm{X}, 0)\bm{A}(\bm{X}, t)\mathrm{d}\bm{X} \\ &= \int_{\partial B_\epsilon^t}\sigma\bm{n}ds(\bm{x}) + \int_{B_\epsilon^t}\bm{f}^{ext}\mathrm{d}\bm{x} \end{aligned} \end{equation}

これは、材料の時間 t の構成におけるすべての B_\epsilon^t \subset \Omega^t に対してである。ここでは、運動量の変化に対する外部体積力(重力など) \bm{f}^{ext} の寄与を加えている。 \bm{f}^{ext} は単位体積あたりの外部体積力を表す。したがって、 B_\epsilon^t 内の運動量の時間変化率は、 B_\epsilon^t に作用する正味の力(コーシー応力場によって表現される表面力と外部体積力の合計)に等しくなる。角運動量保存則のために、 \sigma(\bm{x}, t) は対称でなければならないことも示せる [Bonet and Wood, 2008]。

式(105)の最後の等式は、次のように書くこともできる。

\begin{equation} \int_{B_\epsilon^0}R(\bm{X}, t)J(\bm{x}, t)\bm{A}(\bm{X}, t)\mathrm{d}\bm{X} = \int_{\partial B_\epsilon^t}\sigma\bm{n}\mathrm{d}s(\bm{x}) + \int_{B_\epsilon^t}\bm{f}^{ext}(\bm{x}, t)\mathrm{d}\bm{x} \end{equation}

式(106)の左辺の体積積分をプッシュフォワードすると、両辺が B_\epsilon^t 上の積分になる。

\begin{align} \int_{B_\epsilon^t}\sigma(\bm{x}, t) \bm{a}(\bm{x}, t)\mathrm{d}\bm{x} &= \int_{\partial B_\epsilon^t}\sigma\bm{n}\mathrm{d}s(\bm{x}) + \int_{B_\epsilon^t}\bm{f}^{ext}\mathrm{d}\bm{x} \nonumber\\ &= \int_{B_\epsilon^t}\nabla^{\bm{x}}\cdot\sigma \mathrm{d}\bm{x} + \int_{B_\epsilon^t}\bm{f}^{ext}\mathrm{d}\bm{x} \end{align}

または

\begin{equation} p\bm{a} = \nabla^{\bm{x}}\cdot\sigma+\bm{f}^{ext},\quad\forall\bm{x}\in\Omega^t,t\ge0 \end{equation}

あるいは、式(106)の右辺を次のようにプルバックすることもできる。

\begin{equation} \int_{\partial B_\epsilon^t}\sigma(\bm{x}, t)\bm{n}\mathrm{d}s(\bm{x}) = \int_{\partial B_\epsilon^t}J(\bm{X}, t)\sigma(\phi(\bm{X}, t), t)\bm{F}^{-T}(\bm{X}, t)\bm{N}\mathrm{d}s(\bm{X}) \end{equation}

第一ピオラ-キルヒホッフ応力がコーシー応力と \bm{P} = J\sigma\bm{F}^{-T} として関連していることを思い出すと、次が得られる。

\begin{equation} \int_{\partial B_\epsilon^t}\sigma(\bm{x},t)\mathbf{n}\mathrm{d}s(\bm{x}) = \int_{\partial B_\epsilon^t}\bm{P}(\bm{X},t)\bm{N}\mathrm{d}s(\bm{X}) = \int_{B_\epsilon^t}\nabla^{\bm{x}}\cdot\bm{P}(\bm{X},t)\mathrm{d}\bm{X} \end{equation}

これを式(106)と比較すると、運動量保存則のラグランジュ形式が得られる。

\begin{equation} \int_{B_\epsilon^t}R(\bm{X},0)\bm{A}(\mathbf{X},t)\mathrm{d}\bm{X} = \int_{B_\epsilon^t}\nabla^{\bm{x}}\cdot\bm{P}(\bm{X},t)\mathrm{d}\bm{X} + \int_{B_\epsilon^t}\bm{F}^{ext}J(\bm{X},t)\mathrm{d}\bm{X} \end{equation}

または、

\begin{equation} R(\bm{X},0)\bm{A}(\bm{X},t)=\nabla^{\bm{x}}\cdot\bm{P}(\bm{X},t)+\bm{F}^{ext}(\bm{X},t)J(\bm{X},t),\quad\forall\bm{X}\in\Omega^0,t\geq0 \end{equation}

ここで \bm{F}^{ext} は、単位体積あたりのオイラー体積力 \bm{f}^{ext} のプルバックである。

7-3 弱形式(Weak Form)

MPMは、オイラーグリッド上での応力ベースの力の有限要素法(FEM)による離散化に似ている。したがって、力のつり合いの弱形式(weak form)を使用する。これはラグランジュ的と考えることができる。簡単のため、外力は無視し、運動量保存則から始めることとする。

\begin{equation} R(\bm{X},0)\bm{A}(\bm{X},t) = \nabla^{\bm{x}}\cdot\bm{P}(\bm{X},t),\quad\forall\bm{X},t \end{equation}

または

\begin{equation} R_0A_i=P_{ij,j},\quad\forall\bm{X},t \end{equation}

ここで R_0=R(\bm{X},0) である。任意の関数 \bm{Q}(\cdot,t):\Omega^0\rightarrow\mathbb{R}^d
(テスト関数と呼ばれる)に対して、式(113)との内積を計算し、 \Omega^0 (材料空間)上で積分して弱形式を生成してみる。

\begin{align} \int_{\Omega^0}Q_i(\bm{X},t) R(\bm{X},0)A_i(\bm{X},t)d\bm{X} &= \int_{\Omega^0}Q_i(\bm{X},t) P_{ij,j}(\bm{X},t)\mathrm{d}\bm{X} \\ &= \int_{\Omega^0}((Q_i(\bm{X},t)P_{ij}(\bm{X},t))_{,j} \nonumber \\ &\quad -Q_{i,j}(\bm{X},t) P_{ij}(\bm{X},t))\mathrm{d}\bm{X} \\ &= \int_{\partial \Omega^0}Q_i(\bm{X},t)P_{ij}(\bm{X},t)N_j(\bm{X},t) \mathrm{d}s(\bm{X}) \nonumber \\ &\quad - \int_{\Omega^0}Q_{i,j} (\bm{X},t)P_{ij}(\bm{X},t) \mathrm{d}\bm{X} \end{align}

P_{ij}N_j は境界条件として指定される。単位参照面積あたりの境界力を \bm{T}(\bm{X},t) とし、 T_i = P_{ij}N_j とすると、力のつり合いは、任意の \forall\bm{Q}(\cdot,t):\Omega^0\rightarrow\mathbb{R}^d に対して次を意味すると言える。

\begin{equation} \int_{\Omega^0}Q_i(\bm{X},t)R(\bm{X},0)A_i(\bm{X},t)\mathrm{d}\bm{X} = \int_{\partial\Omega^0}Q_iT_i\mathrm{d}s(\bm{X})-\int_{\Omega^0}Q_{i,j}P_{ij}\mathrm{d}\bm{X} \end{equation}

MPMでは、応力の微分は現在の形状(current configuration)で離散化されるため、応力を含む積分をオイラー表示にプッシュフォワードすることができる。 \bm{Q}(\bm{X},t) = \bm{q}(\phi(\bm{X},t),t) および \bm{q}(\bm{x},t)=\bm{Q}(\phi^{-1}(\bm{x},t),t) となるように \bm{Q} のプッシュフォワードを \bm{q} とすると、次のようになる。

\begin{equation} Q_{i,j}=\frac{\partial Q_i}{\partial X_j}=\frac{\partial q_i}{\partial x_k}\frac{\partial x_k}{\partial X_j}=q_{i,k}F_{kj} \end{equation}

さらに、次を思い出す

\sigma_{ik}=\frac{1}{J}P_{ij}F_{kj}

そして、オイラー形状における単位面積あたりの外力を \bm{t} と定義する( \bm{t}\bm{T} のプッシュフォワード)と、式(118)は次のようになる。

\begin{equation} \int_{\Omega^{\bm{t}}}q_i(\bm{x},t)\rho(\bm{x},t)a_i(\bm{x},t)\mathrm{d}\mathbf{x} = \int_{\partial\Omega^{\bm{t}}}q_it_i\mathrm{d}s(\bm{x}) - \int_{\Omega^{\bm{t}}}q_{i,k}\sigma_{ik}\mathrm{d}\bm{x} \end{equation}

ここでは、セクション5.5で説明した規則を使用して、体積積分と表面積分をプッシュフォワードした。力のつり合いは、上記の式が任意の \bm{q}(\cdot,t):\Omega^t\rightarrow\bm{R}^d に対して成り立つことを意味する。式(120)をオイラー表示における力のつり合いの弱形式と呼ぶ。グリッド上でのMPM離散化は、この方程式に基づく。

8 材料粒子(MATERIAL PARTICLES)

力のつり合いの弱形式(式120)を離散化する前に、まずMPMでシミュレーションされる材料を表す材料粒子(material particles、または材料点 material points、ラグランジュ粒子 Lagrangian particles など)を紹介する。

物質点法(Material Point Method, MPM)は、実際の物質の粒子を追跡するという意味でラグランジュ的であることを思い出す。すなわち、材料粒子の集合 p に対して、質量 (m_p)、速度 (\bm{v}_p)、および位置 (\bm{x}_p) を追跡し続けるのだ。しかし、応力に基づく力はすべてオイラーグリッド上で計算されるため、材料力の効果を取り込むには、材料の状態をオイラー配置(グリッド)に転送する必要がある。その後、これらの効果を材料粒子に戻し、通常のラグランジュ的な方法で粒子を移動させる。このラグランジュ的な性質により、純粋なオイラー法(グリッドベースの流体シミュレーションなど)と比較して、移流(advection)の扱いが非常に簡単になるのだ。

他のPIC/FLIPソルバーと同様に、MPMは背景のオイラーグリッド上で支配方程式を解く。グリッドはスクラッチパッド(計算用の一時的な作業領域)として機能する。各ソルブ(計算ステップ)の後に破棄され、次の時間ステップの開始時に再初期化することができる。実際には、固定された直交格子(fixed Cartesian grid)を使用するのが簡単である。

8-1 オイラー補間関数(Eulerian Interpolating Functions)

MPMの各時間ステップにおいて、粒子はその質量と運動量をグリッドノードに転送している。グリッドでの計算(solve)の後、速度は粒子に移流ステップを実行させるために粒子に戻される。両方の転送には補間関数(interpolation functions)が必要である。有限要素法(FEM)の観点から見ると、オイラーグリッドが本質的な計算メッシュであり、粒子は求積点(quadrature points)として機能する。したがって、補間関数をSPH(Smoothed Particle Hydrodynamics)のように粒子の「カーネル」として考える代わりに、(必須ではないが)補間関数をグリッドノード上で定義されるものとする方がより自然である。

グリッドノード \bm{i} における補間関数を N_{\bm{i}}(\bm{x}) と表記できる。ここで \bm{i} が太字であることに注意する。通常、グリッドノードの多重インデックスだからである。具体的には、2Dでは \bm{i} = (i,j)、3Dでは \bm{i} = (i, k, j) である。 N_{\bm{i}}(\bm{x}) が粒子の位置 \bm{x}_p で評価されるとき、代わりに短い表記 N_{\bm{i}}(\bm{x}_p) = w_{\bm{x}_p} がしばしば使用される。直感的には、各粒子 p とグリッドノード \bm{i} に重み w_{\mathbf{i}p} を関連付けており、これが粒子とノードがどれだけ強く相互作用するかを決定する。もし粒子とグリッドノードが近ければ、重みは大きくなるべきであり、もし粒子とノードが遠く離れていれば、重みは小さくなるべきである。

[Steffen et al., 2008] のように、1次元補間関数のテンソル積(dyadic products)をグリッド基底関数として使用する。

\begin{equation} N_{\mathbf{i}}(\mathbf{x}_p)=N(\frac{1}{h}(x_p-x_{\mathbf{i}}))N(\frac{1}{h}(y_p-y_{\mathbf{i}}))N(\frac{1}{h}(z_p-z_{\mathbf{i}})) \end{equation}

ここで \bm{i}=(i,j,k) はグリッドインデックス、 \bm{x}_p = (x_p, y_p, z_p) は評価位置、 h はグリッド間隔、 \bm{x}_i = (x_i, y_i, z_i) はグリッドノード位置である。よりコンパクトな表記のために、w_{\bm{i}p} = N_{\bm{i}}(\bm{x}_p) および \nabla w_{\bm{i}p} = \nabla N_{\bm{i}}(\bm{x}_p) を使用する。これらをこのように書いている間も、 N_{\bm{i}}\nabla N_{\bm{i}} の両方が依然として任意の \bm{x} の関数であることを心に留めておくべきである。 w で書くときは、これらの関数が \bm{x} = \bm{x}_p で評価されていることを意味する。

カーネル N を選択することは、滑らかさ、計算効率、およびステンシル(影響範囲)の幅に関してトレードオフをもたらす。計算、微分、保存が比較的安価であるため、計算効率の観点からテンソル積スプライン(tensor product splines)を好んで使用する。FLIP流体ソルバーで通常採用される多重線形カーネル(multi-linear kernel)は、これらの選択肢の中で最も単純だが、ここでは適していない。これには二つの理由がある([Steffen et al., 2008]参照)。第一に、 \nabla w_{\bm{i}p} が不連続になり、不連続な力を生成するためであり、第二に、 w_{\bm{i}p}\approx 0 のときでも \nabla w_{\bm{i}p} がゼロから遠い可能性があり、小さな重みを持つグリッドノードに大きな力が加えられる可能性があるためである。MPMは通常、いわゆるセル交差不安定性(cell-crossing instability)を防ぐために、補間関数の C1 連続性(微分が連続であること)を必要とする。実際には、2次Bスプライン(quadratic B splines)または3次Bスプライン(cubic B splines)のいずれかが使用される。2次Bスプラインは、転送ステンシルが小さいため、計算効率が高く、メモリ使用量も節約できる。一方、3次Bスプラインは、より高価だが、より広い範囲をカバーするため、望ましくない芸術的効果でない場合の数値的破壊(numerical fracture)などの数値誤差に対して感度が低くなる。3次カーネルは次のように定義される。

\begin{equation} N(x)= \begin{cases} \frac{1}{2}|x|^3-|x|^2+\frac{2}{3}&0\leq|x|<1 \\ \frac{1}{6}(2-|x|)^3&1\leq|x|<2 \\ 0&2\leq|x| \end{cases} \end{equation}

2次カーネルも有用で

\begin{equation} N(x)= \begin{cases} \frac{3}{4}-|x|^2&0\leq|x|<\frac{1}{2} \\ \frac{1}{2}(\frac{3}{2}-|x|)^2&\frac{1}{2}\leq|x|<\frac{3}{2} \\ 0&\frac{3}{2}\leq|x| \end{cases} \end{equation}

2次および3次カーネルを図3に示す 。理論的には不安定ですが、線形補間関数も実際には特定の用途で機能する可能性がある 。それは、

\begin{equation} N(x) = \begin{cases} 1 - |x| & 0 \le |x| < 1 \\ 0 & 1 \le |x| \end{cases} \end{equation}

勾配の計算は、同様に1次元関数を微分することによって行われる。

\nabla N_{\bm{i}}(\bm{x}_p) = \begin{pmatrix} \frac{1}{h} N' \! \left( \frac{1}{h} (x_p - x_{\bm{i}}) \right) N \! \left( \frac{1}{h} (y_p - y_{\bm{i}}) \right) N \! \left( \frac{1}{h} (z_p - z_{\bm{i}}) \right) \\ % 1行目 N \! \left( \frac{1}{h} (x_p - x_{\bm{i}}) \right) \frac{1}{h} N' \! \left( \frac{1}{h} (y_p - y_{\bm{i}}) \right) N \! \left( \frac{1}{h} (z_p - z_{\bm{i}}) \right) \\ % 2行目 N \! \left( \frac{1}{h} (x_p - x_{\bm{i}}) \right) N \! \left( \frac{1}{h} (y_p - y_{\bm{i}}) \right) \frac{1}{h} N' \! \left( \frac{1}{h} (z_p - z_{\bm{i}}) \right) % 3行目 \end{pmatrix}

ここで N'(x)N(x) の導関数である。

8-2 オイラー/ラグランジュ質量(Eulerian/Lagrangian Mass)

我々の材料を有限な点の集合として表現する際、各点に全材料の一部であるサブセット B_{\Delta x,p}^0\subset\Omega^0 を割り当てることができる。このようにして、その粒子の質量を次のように定義できる。

\begin{equation} m_p^n=\int_{B_{\Delta x,p}^{t^n}}\rho(\bm{x},t^n)\mathrm{d}\bm{x} \end{equation}

この慣例を用いると、質量と運動量(そして次に速度)をオイラーグリッドのノードに転送するための、以下の保存的なプロセスを定義でき、次のように定義する。

\begin{equation} m_{\bm{i}}=\sum_pm_pN_{\bm{i}}(\bm{x}_p) \end{equation}

これをオイラーグリッドノード \bm{i} の質量とする。この慣例を用いると、次のようになる。

\begin{equation} \sum_{\bm{i}}m_{\bm{i}} = \sum_{\bm{i}}\sum_pm_pN_{\bm{i}}(\bm{x}_p) = \sum_pm_p\sum_{\bm{i}}N_{\bm{i}}(\bm{x}_p) = \sum_pm_p \end{equation}

これは、 N_{\bm{i}} に対する「1の分割(partition of unity)」の仮定によるのもである。

8-3 オイラー/ラグランジュ運動量(Eulerian/Lagrangian Momentum)

同様に、粒子の運動量 m_p\bm{v}_p を転送できる(注意:最終的にはセクション10.1で説明されるAPICと呼ばれる異なる転送方法を使用する)。

\begin{equation} (m\bm{v})_{\bm{i}} = \sum_pm_p\bm{v}_pN_{\bm{i}}(\bm{x}_p) \end{equation}

そして、次のことを示すことができる。

\begin{equation} \sum_{\bm{i}}(m\bm{v})_{\bm{i}} = \sum_pm_p\bm{v}_p \end{equation}

これは、同じ「1の分割(partition of unity)」のロジックによるものである。

オイラー速度 \bm{v}_{\bm{i}} は次のように定義される。

\begin{equation} \bm{v}_{\bm{i}} = \frac{(m\bm{v})_{\bm{i}}}{m_{\bm{i}}} \end{equation}

この場合、繰り返された添字 \bm{i} は合計を意味しないことに注意する。

8-4 オイラーからラグランジュへの転送(Eulerian to Lagrangian Transfer)

オイラー変数(グリッド上の変数)からラグランジュ変数(粒子上の変数)への転送も同様である。しかし、ラグランジュ粒子の質量は決して変化しないため、グリッドから粒子へ質量を転送する必要はないのである。速度は単純に次のように補間される。

\begin{equation} \bm{v}_p=\sum_{\bm{i}}\bm{v}_{\bm{i}}N_{\bm{i}}(\bm{x}_p) \end{equation}

これが運動量を保存することは、次のように容易に検証できる。

\begin{equation} \sum_{p}m_p\bm{v}_p=\sum_pm_p\sum_{\bm{i}}\bm{v}_{\bm{i}}N_{\bm{i}}(\bm{x}_p)=\sum_{\bm{i}}\bm{v}_{\bm{i}}\sum_pm_pN_{\bm{i}}(\bm{x}_p)=\sum_{\bm{i}}m_{\bm{i}}\bm{v}_{\bm{i}} \end{equation}

9 離散化(DISCRETIZATION)

9-1 離散時間(Discrete Time)

まず、時刻 t^n にいると仮定することから始める。力のつり合い方程式の弱形式(式(118)および(120))は、次を意味することを思い出す。

\begin{equation} \int_{\Omega^0}Q_iR_0A_id\bm{X}=\int_{\partial\Omega^{t^n}}q_it_i\mathrm{d}s(\bm{x})-\int_{\Omega^{t^n}}q_{i,k}\sigma_{ik}\mathrm{d}\bm{x} \end{equation}

すべての \bm{q}(\bm{x},t^n)(または \bm{Q}(\bm{X}, t^n))に対して。まず、ラグランジュ加速度 A_i(\bm{X},t^n)\frac{1}{\Delta t}\left(V_i^{n+1} - V_i^n\right) で置き換えることから始める。次に、左辺を \Omega^0 (材料空間)から \Omega^{t^n}(現在の空間)へプッシュフォワードして、次を得ることができる。

\begin{align} &\frac{1}{\Delta t}\int_{\Omega^{t^n}}q_i(\bm{x}, t^n)\rho(\bm{x},t^n)\left(v_i^{n+1}(\bm{x})-v_i^n(\bm{x})\right)\mathrm{d}\bm{x} \nonumber \\ % 2行目: = と右辺の最初の項 &= \int_{\partial\Omega^{t^n}}q_i(\bm{x},t^n)t_i(\bm{x},t^n)\mathrm{d}s(\bm{x})-\int_{\Omega^{t^n}}q_{i,k}(\bm{x}, t^n)\sigma_{ik}(\bm{x},t^n)\mathrm{d}\bm{x} \end{align}

この表記では、 \bm{v}^n:\Omega^{t^n}\rightarrow\mathbb{R}^d,\bm{v}^{n+1}:\Omega^{t^n}\rightarrow\mathbb{R}^d(両方とも \bm{x}\in\Omega^{t^n} に対して定義される)であることに注意する。したがって、 v_i^{n+1}(\bm{x}) = V_i(\phi^{-1}(\bm{x},t^n),t^{n+1}) であり、 v_i^n(\bm{x}) = V_i(\phi^{-1}(\bm{x},t^n),t^n) である。

式(134)における添字 ik はすべて、次元の成分インデックスであることに留意すべきである。すなわち、3Dでは i = 1, 2, 3、2Dでは i = 1, 2 である。空間を離散化する際には、グリッドノードを示すために i,j,k のような添字を使いたくなるでしょう。混乱を避けるために、これ以降、成分インデックスを示すためにはギリシャ文字、例えば \alpha,\beta,\gamma を必要とすることができる。これにより、式(134)は次のように書き換えられる。

\begin{align} &\frac{1}{\Delta t}\int_{\Omega^{t^n}}q_\alpha(\bm{x},t^n)\rho(\bm{x},t^n)\left(v_\alpha^{n+1}(\bm{x})-v_\alpha^n(\bm{x})\right)\mathrm{d}\bm{x} \nonumber \\ % 2行目: = と右辺の最初の項 &= \int_{\partial\Omega^{t^n}}q_\alpha(\bm{x},t^n)t_\alpha(\bm{x},t^n)\mathrm{d}s(\bm{x})-\int_{\Omega^{t^n}}q_{\alpha,\beta}(\bm{x}, t^n)\sigma_{\alpha\beta}(\bm{x},t^n)\mathrm{d}\bm{x} \end{align}

ここで、 \alpha,\beta = 1,2,...d であり、 d は問題の次元(2または3)である。

そして、 q_{\bm{i}\alpha} は、ノード \bm{i} に格納されているベクトル量 \bm{q}\alpha 成分を意味する。また、 q_\alpha(\bm{x},t) は、(連続的な)場 q_\alpha(\bm{x},t)\alpha 成分を意味する。

9-2 離散空間(Discrete Space)

空間項の有限要素法(FEM)スタイルの離散化は、 q_\alpha,v_\alpha^n,v_\alpha^{n+1} をグリッドベースの補間表現で置き換えることによって行うことができる。

\begin{gather} q_\alpha(\bm{x},t^n) = \sum_{\bm{i}}q_{\bm{i}\alpha}(t^n)N_{\bm{i}}(\bm{x}) \\ v_\alpha^n(\bm{x}) = \sum_{j}v_j{\alpha}^nN_{j}(\bm{x}) \\ v_\alpha^{n+1}(\bm{x}) = \sum_jv_{j\alpha}^{n+1}N_j(\bm{x}) \end{gather}

または、短く書くと(繰り返された添字は合計を意味する)

\begin{equation} q_\alpha^n = q_{\bm{i}\alpha}^nN_{\bm{i}},\quad v_\alpha^n=v_{\bm{j}\alpha}^nN_{\bm{j}},\quad v_\alpha^{n+1} = v_{\bm{j}\alpha}^{n+1}N_{\bm{j}} \end{equation}

力のつり合い(式135)は、次のように見なすことができる。すべての q_{\bm{i}\alpha} に対して

\begin{align} % 1行目: 左辺の第1項 &\frac{1}{\Delta t}\int_{\Omega^{t^n}}q_{\bm{i}\alpha}^nN_{\bm{i}}(\bm{x})\rho(\bm{x}, t^n)v_{\bm{j}\alpha}^{n+1} N_{\bm{j}}(\bm{x})\,\mathrm{d}\bm{x} \nonumber\\ &\qquad\qquad\qquad-\frac{1}{\Delta t}\int_{\Omega^{t^n}}q_{\bm{i}\alpha}^n N_{\bm{i}}(\bm{x})\rho(\bm{x}, t^n)v_{\bm{j}\alpha}^nN_{\bm{j}}(\bm{x})\,\mathrm{d}\bm{x}\nonumber \\ % 3行目: = と右辺の第1項 &= \int_{\partial\Omega^{t^n}}q_{\bm{i}\alpha}^nN_{\bm{i}}(\bm{x})t_\alpha(\bm{x},t^n)\,\mathrm{d}s(\bm{x})\nonumber \\ &\qquad\qquad\qquad- \int_{\Omega^{t^n}} q_{\bm{i}\alpha}^n N_{\bm{i},\beta}(\bm{x})\sigma_{\alpha\beta}(\bm{x},t^n)\, \mathrm{d}\bm{x} \end{align}

である。これは次のように表現することもできる。

\begin{align} &q_{i\alpha}^n\delta_{\alpha\beta}\frac{m_{\bm{i}\bm{j}}^n}{\Delta t}v_{\bm{j}\beta}^{n+1} - q_{\bm{i}\alpha}^n\delta_{\alpha\beta}\frac{m_{\bm{i}\bm{j}}^n}{\Delta t}v_{\bm{j}\beta}^n \nonumber \\ &= \int_{\partial\Omega^{t^n}}q_{\bm{i}\alpha}^nN_{\bm{i}}t_{\alpha}\mathrm{d}s(\bm{x}) - \int_{\Omega^{t^n}}q_{\bm{i}\alpha}^nN_{\bm{i},\beta}\sigma_{\alpha\beta}\mathrm{d}\bm{x} \end{align}

ここで、

\begin{equation} m_{\bm{i}\bm{j}}^n = \int_{\Omega^{t^n}}N_{\bm{i}}(\bm{x})\rho(\bm{x},t^n)N_{\bm{j}}(\bm{x})\mathrm{d}\bm{x} \end{equation}

は質量行列(mass matrix)を定義している。これを材料空間にプルバックすると、次が得られる。

\begin{equation} m_{\bm{i}\bm{j}}^n = \int_{\Omega^0}N_{\bm{i}}(\bm{x}(\bm{X}))R(\bm{X},0)N_{\bm{j}}(\bm{x}(\bm{X}))d\bm{X} \approx \sum_pm_pN_{\bm{i}}(\bm{x}_p)N_{\bm{j}}(\bm{x}_p) \end{equation}

これは、 B_{ip}=\sqrt{m_p}N_{ip} とおくと m_{ij} = \Sigma_pB_{ip}B_{jp} (\bm{B}\bm{B}^T の形) と書けるため、対称正定値(または半正定値)である。。したがって、任意の \bm{z} に対して \bm{z}^T\bm{M}\bm{z}\ge0 である。実際には、この質量行列は特異(singular、逆行列が存在しない)になる可能性があるため、使用できない。以下で、粒子-グリッド転送と整合性を保ちながら、質量行列を対角かつ正定値な行列で近似する「集中質量(mass lumping)」戦略を見ていく。

[](data:image/svg+xml;utf8,<svg xmlns="http://www.w3.org/2000/svg" width="400em" height="1.08em" viewBox="0 0 400000 1080" preserveAspectRatio="xMinYMin slice"><path d="M95,702
c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,-10,-9.5,-14
c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54
c44.2,-33.3,65.8,-50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10
s173,378,173,378c0.7,0,35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429
c69,-144,104.5,-217.7,106.5,-221
l0 -0
c5.3,-9.3,12,-14,20,-14
H400000v40H845.2724
s-225.272,467,-225.272,467s-235,486,-235,486c-2.7,4.7,-9,7,-19,7
c-6,0,-10,-1,-12,-3s-194,-422,-194,-422s-65,47,-65,47z
M834 80h400000v40h-400000z"></path></svg>)

式(141)は、 q_{\bm{i}\alpha}^n のすべての選択に対して真でなければならない。そこで、次のように選択する

q_{\bm{i}}^n = \begin{cases} 1,\quad\alpha = \hat{\alpha}~\text{and}~\bm{i}=\hat{\bm{i}} \\ 0,\quad\text{otherwise} \end{cases}

すると次のようになる。

\begin{equation} \sum_{\bm{j}}\frac{m_{\hat{\bm{i}\bm{j}}}}{\Delta t}\left(v_{\bm{j}\hat{\alpha}}^{n+1}-v_{\bm{j}\hat{\alpha}}^{n}\right)=\int_{\partial\Omega^{t^n}}N_{\hat{\bm{i}}}t_{\hat{\alpha}}\mathrm{d}s(\bm{x}) - \int_{\Omega^{t^n}}N_{\hat{\bm{i}},\beta}\sigma_{\hat{\alpha}\beta}\mathrm{d}\bm{x} \end{equation}

これは、グリッド上の \hat{\bm{i}},\hat{\alpha} 自由度に対する離散的な力のつり合い方程式と見なすことができる。次に、これがオイラー運動量の陽的な更新規則として使用でき、右辺が \hat{\bm{i}}^{th} 番目のオイラーグリッドノードに作用する力の \hat{\alpha}^{th} 成分と考えられることを示す。

先述のように、集中質量(mass lumping)による単純化を使用すると便利なことが多い。(ただし精度は低下する)。これは、 m_{\bm{i}\bm{j}} の行を対応する行和で置き換え、それによって対角行列にすることによって行われる。行和(行 \bm{i} に対して \hat{m}_{\bm{i}} と呼ぶ)は次のようになる。

\begin{align} % 1行目: 左辺と最初の右辺 \hat{m}_{\bm{i}} &= \sum_{\bm{j}} \int_{\Omega^{t^n}} N_{\bm{i}}(\bm{x}) \rho(\bm{x}, t^n) N_{\bm{j}}(\bm{x}) \, \mathrm{d}\bm{x} \nonumber \\ % 2行目: 2番目の右辺 (= で揃える) &= \int_{\Omega^{t^n}} N_{\bm{i}}(\bm{x}) \rho(\bm{x}, t^n) \sum_{\bm{j}} N_{\bm{j}}(\bm{x}) \, \mathrm{d}\bm{x} \nonumber \\ % 3行目: 最後の右辺 (= で揃える) と式番号 &= \int_{\Omega^{t^n}} N_{\bm{i}}(\bm{x}) \rho(\bm{x}, t^n) \, \mathrm{d}\bm{x} \end{align}

次の近似(式143で行ったことと同様)も使用できる。

\begin{equation} \int_{\Omega^{t^n}} N_{\bm{i}}(\bm{x}) \rho(\bm{x}, t^n) \, \mathrm{d}\bm{x} = \int_{\Omega^0} N_{\bm{i}}(\bm{x}(\bm{X})) R(\bm{X}, 0) \, \mathrm{d}\bm{X} \approx \sum_p N_{\bm{i}}(\bm{x}_p) m_p \end{equation}

言い換えれば、 \hat{m}_{\bm{i}} は自然に m_{\bm{i}}(式(126)で定義されたグリッドノード質量)として近似できる。

m_{\bm{i}}\bm{v}_{\bm{i}}^n を運動量 (m\bm{v})_{\bm{i}}^n として書き直すと、離散化を次のように要約できる。

\begin{align} % 左辺 (分数) \frac{(mv)_{i\alpha}^{n+1} - (mv)_{i\alpha}^n}{\Delta t} % 右辺 (積分項の差) = \int_{\partial \Omega^{t^n}} &N_i(x) t_\alpha(x, t^n) \, \mathrm{d}s(x)\nonumber \\& - \int_{\Omega^{t^n}} N_{i,\beta}(x) \sigma_{\alpha\beta}(x, t^n) \, \mathrm{d}x \end{align}

言い換えれば、左辺は運動量の変化であるため、右辺はおおよそ力である。

ここでは、各ラグランジュ粒子 x_p^n における応力の推定値 \sigma_p^n\approx\sigma(x_p^n,t^n) があると仮定した。したがって、

\begin{equation} \int_{\Omega^{t^n}}N_{\bm{i},\beta}(\bm{x})\sigma_{\alpha\beta}(\bm{x},t^n)\mathrm{d}\bm{x}\approx\sum_p{\sigma_{p}^n}_{\alpha\beta}N_{\bm{i},\beta}(\bm{x}_p^n)V_p^n \end{equation}

ここで V_p^nB_{\Delta x,p}^{t^n} の体積、または粒子 p が時刻 t^n で占める体積である。

9-3 体積の推定(Estimating the Volume)

V_p^n(粒子 p の時刻 t における体積)を推定するには、二つの方法がある。第一に、次を思い出す。

\begin{equation} m_p\approx R(\bm{X}_p,0)V_p^0\approx\rho(\bm{x}_p^n,t^n)V_p^n \end{equation}

次に、 m_{\bm{i}}^n から \rho(\bm{x}_p^n,t^n) を近似することができる。

\begin{align} \rho_{\bm{i}}^n&=\frac{m_{\bm{i}}^n}{\Delta x^d} \\ \rho(\bm{x}_p^n,t^n)&\approx\sum_{\bm{i}}\rho_{\bm{i}}^nN_{\bm{i}}(\bm{x}_p^n) \end{align}

従って、

\begin{equation} V_p^n\approx\frac{m_p}{\Sigma_{\bm{i}}\frac{m_{\bm{i}}^n}{\Delta x^d}N_{\bm{i}}(\bm{x}_p^n)} = \frac{m_p\Delta x^d}{\Sigma_{\bm{i}}m_{\bm{i}}^nN_{\bm{i}}(\bm{x}_p^n)} \end{equation}

別のアプローチは次の通りである。各ラグランジュ粒子における変形の近似 \bm{F}_p^n\approx\bm{F}(\bm{X}_p,t^n) がある場合(これは通常、弾性材料で必要である)、 J_p^n = \mathrm{det}(\bm{F}_p^n) を用いて V_p^n を近似できる。具体的には、初期体積 V_p^0 = \int_{B_{\Delta x,p}^0} を知っていると仮定すると、

\begin{equation} V_p^n\approx V_p^0J_p^n \end{equation}

第一ピオラ-キルヒホッフ応力の式 \sigma = \frac{1}{J}\bm{P}\bm{F}^T を使用する場合、式(147)および(148)を、 \sigma の代わりに \bm{P} を用いてさらに書き換えることができる。

\begin{align} \sum_p{\sigma_{p}^n}_{\alpha\beta}N_{\bm{i},\beta}(\bm{x}_p^n)V_p^n &= \sum_{p}\frac{1}{J_p^n}{P_p^n}_{\alpha\gamma}{F_p^n}_{\beta\gamma}N_{\bm{i},\beta}(\bm{x}_p^n)V_p^0J_p^n \nonumber\\ &= \sum_p{P_p^n}_{\alpha\gamma}{F_p^n}_{\beta\gamma}N_{\bm{i},\beta}(\bm{x}_p^n)V_p^0 \end{align}

この式は、ほとんどの構成モデルが \bm{P} で表現されるため、実際にはより有用である。さて、

\begin{align} \frac{\left((mv)_{\bm{i}\alpha}^{n+1} - (mv)_{\bm{i}\alpha}^n\right)}{\Delta t} = \int_{\partial\Omega^{t^n}}&N_{\bm{i}}(\bm{x})t_\alpha(\bm{x},t^n)\mathrm{d}s(\bm{x}) \nonumber\\ &- \sum_p{P_p^n}_{\alpha\gamma}{F_p^n}_{\beta\gamma}N_{\bm{i},\beta}(\bm{x}_p^n)V_p^0 \end{align}

次に、各ラグランジュ点( x_p^n )における変形勾配( F_p^n )をどのように近似するかについて議論する。

9-4 変形勾配の発展(Deformation Gradient Evolution)

時間ステップスキームの議論を続けるためには、構成モデル(constitutive model)を考慮する必要がある。これまでの結果は材料の種類に依存しないものであった。我々は一般的に、応力が主に変形勾配によって表現される材料の形状変化に依存するモデルを考える(これにはセクション6で議論したすべてのモデルが含まれる)。次に、これをMPMの文脈でどのように計算できるかについて議論をする。変形勾配は、オイラー法(またはメッシュがないラグランジュ法、例: 粒子法)では計算がより困難である。しかし、方程式(セクション5.4の結果を思い出してください)

\begin{equation} \frac{\partial}{\partial t}\bm{F}(\bm{X},t) = \frac{\partial \bm{v}}{\partial \bm{x}}(\phi(\bm{X},t),t)\bm{F}(\bm{X},t) \end{equation}

を用いて、オイラーグリッド上で \frac{\partial \bm{v}}{\partial \bm{x}} (速度勾配)を離散化することにより、各材料粒子状の変形勾配を更新(または発展)させることができる。さて、同様のルールを離散的に再定式化する。時刻 t^{n+1} における速度( \Omega^{t^n} 上で定義される)が \bm{x} 飲みの関数として得られていると仮定する。具体的には、 \bm{v}^{n+1}(\bm{x})=\bm{V}(\phi^{-1}(\bm{x},t^n),t^{n+1}) と定義される関数 \bm{v}^{n+1}:\Omega^{t^n}\rightarrow\mathbb{R}^d を考える。もちろん \bm{v}^{n+1}(\phi(\bm{X},t^n)) = \bm{V}(\bm{X},t^{n+1}) である。(セクション9.1で同じ定義を使用した)これを用いると、次のようになる。

\begin{equation} \frac{\partial}{\partial t}\bm{F}(\bm{X}, t^{n+1}) = \frac{\partial \bm{V}}{\partial \bm{X}}(\bm{X},t^{n+1}) = \frac{\partial \bm{v}^{n+1}}{\partial \bm{x}}(\phi(\bm{X},t^n))\bm{F}(\bm{X},t^n) \end{equation}

これは便利である。なぜなら、さらに次の近似を用いると、

\begin{equation} \frac{\partial}{\partial t}\bm{F}(\bm{X}_p,t^{n+1})\approx\frac{\bm{F}_p^{n+1}-\bm{F}_p^n}{\Delta t} \end{equation}

次のようになる。

\begin{equation} \bm{F}_p^{n+1} = \bm{F}_p^n+\Delta t\frac{\partial \bm{v}^{n+1}}{\partial \bm{x}}(\bm{x}_p^n)\bm{F}_p^n = \left(\bm{I}+\Delta{t}\frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{x}_p^n)\right) \end{equation}

そして、もちろん \bm{v}^{n+1} に対してグリッドベースの補間式、すなわち

\begin{align} \bm{v}^{n+1}(\bm{x}) &= \sum_{\bm{i}}\bm{v}_{\bm{i}}^{n+1}N_{\bm{i}}(\bm{x}) \\ \frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{x}) &= \sum_{\bm{i}}\bm{v}_{\bm{i}}^{n+1}\left(\frac{\partial N_{\bm{i}}}{\partial\bm{x}}(\bm{x})\right)^T \end{align}

を用いると次のようになる。

\begin{equation} \bm{F}_p^{n+1} = \left(\bm{I}+\Delta{t}\sum_{\bm{i}}\bm{v}_{\bm{i}}^{n+1}\left(\frac{\partial N_{\bm{i}}}{\partial\bm{x}}(\bm{x}_p^n)\right)^T\right)\bm{F}_p^n \end{equation}

これが、 \bm{v}_{\bm{i}}^{n+1}\bm{F}_p^n が与えられた時の \bm{F}_p^{n+1} の更新ルールとなる。

9-5 エネルギー勾配としての力(Forces as Energy Gradient)

弾性応答は、弾性ポテンシャルエネルギーから生じることが示せる。これは連続的な場合と離散的な場合の両方で真である。離散的な場合について、 \hat{\bm{x}}_i = \bm{x}_i であるときのオイラーグリッドのノードを \hat{\bm{x}}_i と定義することを考える。言い換えれば、一時的に、オイラーグリッドのノードが、仮想的な移動後のノード位置を示す変数 \hat{\bm{x}}_i で定義される変数で動くことを許容したと考える。すると、これら \hat{\bm{x}}_i を一時的に新しいラグランジュ自由度と考えることができる。これは、ラグランジュ粒子として \bm{X}_p を持つことから切り替えて、 \bm{X}_{\bm{i}} = \phi^{-1}(\bm{x}_{\bm{i}},t^n) を新しいラグランジュ粒子とし、これによって \bm{x}_{\bm{i}}^{n+1} = \phi(\bm{X}_{\bm{i}},t^{n+1}) = \hat{\bm{x}}_{\bm{i}} が材料の新しい構成を定義すると言うようなものである。この考え方を念頭に置くと、セクション9.2で導出された力が、ある種の弾性材料の場合には離散的なポテンシャルエネルギーに関連していることを示すことができる。第一ピオラ-キルヒホッフ応力が弾性ポテンシャルエネルギー密度 \psi と次のように関連している超弾性材料に対して、

\begin{equation} \bm{P}(\bm{F}) = \frac{\partial \psi}{\partial \bm{F}}(\bm{F}) \end{equation}

セクション6で見たように、全ポテンシャルエネルギー(移動したノード位置 \hat{\bm{x}} の関数として)を次のように定義できる。

\begin{equation} e(\hat{\bm{x}}) = \sum_p \psi(\bm{F}_p(\hat{\bm{x}}))V_p^0 \end{equation}

ここで \hat{\bm{x}} はすべてのオイラー \hat{\bm{x}}_{\bm{i}} の完全な(組み立てられた)ベクトルである。ここでは、変形勾配を \hat{\bm{x}} の関数と考える。これは、一時的に材料の運動を \bm{X}_{\bm{i}} = \phi^{-1}(\bm{x}_{\bm{i}},t^n)\hat{\bm{x}}_{\bm{i}} で定義されると考えているからである。実際には、これは運動を \bm{v}_{\bm{i}}^{n+1} で定義することと同等の概念である。具体的には、次のように定義されると考えることができる。

\begin{equation} \bm{v}_{\bm{i}}^{n+1} = \frac{\hat{\bm{x}}_{\bm{i}} - \bm{x}_{\bm{i}}}{\Delta t}\quad\text{or}\quad\hat{\bm{x}}_{\bm{i}} = \bm{x}_{\bm{i}}+\Delta t\bm{v}_{\bm{i}}^{n+1} \end{equation}

これは、粒子ごとの変形勾配の式につながる。

\begin{equation} {F_p}_{\beta\gamma}(\hat{\bm{x}}) = {F_p^n}_{\beta\gamma}+\Delta t\sum_{\bm{j}}\left(\frac{\hat{x}_{\bm{j}\beta} - x_{\bm{j}\beta}}{\Delta t}\right)\sum_{\tau}N_{\bm{j},\tau}(\bm{x}_p^n){F_p^n}_{\tau\gamma} \end{equation}

エネルギーを \hat{\bm{x}}_{\bm{i}\alpha}\bm{X}_{\bm{i}} [訳注: \hat{\bm{x}}_{\bm{i}}]の \alpha 成分)について微分すると、次のようになる。

\begin{equation} \frac{\partial e}{\partial\hat{x}_{\bm{i}\alpha}}(\hat{\bm{x}}) = \sum_p\sum_{\beta,\gamma} P_{\beta\gamma}(\bm{F}_p(\hat{\bm{x}}))\frac{\partial {F_p}_{\beta\gamma}}{\partial\hat{x}_{\bm{i}\alpha}}(\hat{\bm{x}})V_p^0 \end{equation}

上記の変形勾配の式から、次がわかる。

\begin{equation} \frac{\partial {F_p}_{\beta\gamma}}{\partial\hat{x}_{\bm{i}\alpha}}(\hat{\bm{x}}) = \delta_{\alpha\beta}\sum_{\tau}N_{\bm{i},\tau}(\bm{x}_p^n){F_p^n}_{\tau\gamma} \end{equation}

そして、これを式(167)に代入すると、次が得られる。

\begin{equation} \frac{\partial e}{\partial \hat{x}_{\bm{i}\alpha}}(\hat{\bm{x}}) = \sum_pP_{\alpha\gamma}(\bm{F}_p(\hat{\bm{x}})){F_p^n}_{\tau\gamma}N_{\bm{i},\tau}(\bm{x}_p^n)V_p^0 \end{equation}

したがって、セクション9.3 [訳注: 9.2] の結果と比較すると、オイラーグリッドノード \bm{i} に作用する力は次のように見ることができる。

\begin{equation} -\frac{\partial e}{\partial\hat{x}_{\bm{i}\alpha}}\left(\hat{\bm{x}} = \begin{pmatrix} \bm{x}_0 \\ \bm{x}_1 \\ \vdots \end{pmatrix}\right) \end{equation}

言い換えれば、運動量の更新を次のように書くことができる。

\begin{align} (mv)_{\bm{i}\alpha}^{n+1} = (mv)_{\bm{i}\alpha}^n &- \Delta t\frac{\partial e}{\partial \hat{x}_{\bm{i}\alpha}}\left(\hat{\bm{x}} = \begin{pmatrix} \bm{x}_0 \\ \bm{x}_1 \\ \vdots \end{pmatrix} \right) \nonumber\\&+ \Delta t\int_{\partial\Omega^{t^n}}N_{\bm{i}}(\bm{x})t_\alpha(\bm{x},t^n)\mathrm{ds}(\bm{x}) \end{align}

10 陽的時間積分(EXPLICIT TIME INTEGRATION)

このセクションでは、シンプレクティック・オイラー(Symplectic Euler)時間積分を用いたMPMの最も簡単な実装について説明する。このセクションのほとんどの数式は、前のセッションで導出されたものである。ここでは、MPMを実装するためのより簡潔で実践者にやさしいリファレンスを提供するために、それらを再喝する。まずいくつかの個別のモジュールについて議論し、次に完全なMPMアルゴリズムを述べる。まずいくつかの個別のモジュールについて議論し、次に完全なMPMアルゴリズムを述べる。

10-1 APIC転送(APIC Transfers)

アフィンPIC(Affine Particle-In-Cell, APIC)法 [Jiang et al., 2015] が、その良好な数値特性のため、粒子-グリッド間転送に推奨される。ここでは表記の簡単のため、時間の上付き添字は無視する。

粒子からグリッドへの転送は次のようになる。

\begin{align} m_{\bm{i}} &= \sum_pw_{ip}m_p \\ m_{\bm{i}}\bm{v}_{\bm{i}} &= \sum_pw_{ip}m_p(\bm{v}_p+\bm{B}_p(\bm{D}_p)^{-1}(\bm{x}_{\bm{i}} - \bm{x}_p)) \end{align}

ここで、 \bm{B}_p は各粒子に(質量、位置、速度と同様に)格納される行列量であり、 \bm{D}_p は次で与えられる。

\begin{equation} \bm{D}_p = \sum_{\bm{i}}w_{ip}(\bm{x}_{\bm{i}}-\bm{x}_p)(\bm{x}_{\bm{i}}-\bm{x}_p)^T \end{equation}

そして、転送中にアフィン運動(平行移動、回転、スケール、せん断を含む線形変換)を保存するように導出される。対応するグリッドから粒子への転送は次のようになる。

\begin{align} \bm{v}_p &= \sum_{\bm{i}}w_{ip}\bm{v}_{\bm{i}} \\ \bm{B}_p &= \sum_{\bm{i}}w_{ip}\bm{v}_{\bm{i}}(\bm{x}_{\bm{i}}-\bm{x}_p)^T \end{align}

都合の良いことに、 \bm{D}_p はMPMで一般的に使用される2次補間ステンシル( \bm{D}_p = \frac{1}{4}\Delta x^2\bm{I} )および3次補間ステンシル( \bm{D}_p = \frac{1}{3}\Delta x^2\bm{I} )の場合、驚くほど単純な形をとる。これらの補間ステンシルでは、 (\bm{D}_p)^{-1} を掛けることは定数倍のスケーリングに相当することに注意する。

もし、(推奨されないが)線形スプラインが補間関数として使用される場合、 \bm{D}_p が特異になるのを防ぐために、 APICの粒子からグリッドへの代替定式化が次のように与えられる。

\begin{align} m_{\bm{i}} = \sum_pw_{ip}m_p \\ m_{\bm{i}}\bm{v}_{\bm{i}} = \sum_p w_{ip}m_p(\bm{v}_p+\bm{C}_p(\bm{x}_{\bm{i}}-\bm{x}_p)) \end{align}

そして、グリッドから粒子へは次のようになる。

\begin{align} \bm{v}_p &= \sum_{\bm{i}}w_{ip}\bm{v}_{\bm{i}} \\ \bm{C}_p &= \sum_{\bm{i}}\bm{v}_{\bm{i}}\left(\frac{\partial N_{\bm{i}}}{\partial \bm{x}}(\bm{x}_p)\right)^T \end{align}

この場合、 \bm{C}_p は単に \bm{x}_p で評価されたオイラー速度勾配であり、 \bm{F}_p を発展させるために既に計算したものであることに注意する。

10-2 変形勾配の更新(Deformation Gradient Update)

各時間ステップの開始時に置いて、グリッドノードの位置は常に変形していない正則なグリッド上にある。それらを \bm{x}_{\bm{i}}^n と呼ぶこととする。常に同じグリッドを使用する場合は上付き添え字 n は不要だが、ここでは表記の慣例として残しておく。各粒子の変形勾配は、グリッドの運動に従って更新することができる。グリッド速度場 \bm{v}_{\bm{i}}^{n+1} を仮定すると、 \bm{F} の更新ルールは次のように示せる。

\begin{equation} \bm{F}_p^{n+1} = \left(\bm{I} + \Delta t\sum_{\bm{i}}\bm{v}_{\bm{i}}^{n+1}(\nabla w_{\bm{i}p}^n)^T\right)\bm{F}_p^n \end{equation}

10-3 状態の更新(State Update)

グリッド上でのシンプレクティック・オイラー時間積分は、次を意味する。

\begin{align} \bm{x}_{\bm{i}}^{n+1} &= \bm{x}_{\bm{i}}^n + \Delta t\bm{v}_{\bm{i}}^{n+1} \\ \bm{v}_{\bm{i}}^{n+1} &= \bm{v}_{\bm{i}^n}+\Delta t \mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^n)/m_{\bm{i}} \end{align}

グリッドの運動についての一つの考え方は、

\begin{equation} \bm{x} = \bm{x}(\bm{v}) \end{equation}

すなわち、グリッドノードの位置は単純に式(182)を介してグリッドノード速度の関数である、というものである。従って、

\begin{equation} \mathbf{f}_{\bm{i}}^n = \mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^n) = \mathbf{f}_{\bm{i}}(\bm{x}(0)) \end{equation}

すなわち、陽的な力(explicit force)とは、単純にグリッドの運動がゼロであると仮定した場合の力である。式(182)を用いると、式(181)を次のように書き換えることができる。

\begin{equation} \bm{F}_p(\bm{x}) = \left(\bm{I} + \sum_{\bm{i}}(\bm{x}_{\bm{i}} - \bm{x}_{\bm{i}}^n)(\nabla w_{\bm{i}p}^n)^T\right)\bm{F}_p^n \end{equation}

これにより、 \bm{F}_p はグリッドノード位置 \bm{x} の関数となる。

10-4 力(Forces)

MPMにおける力はグリッドノード上で定義される。以前に見たように、超弾性固体の場合、力は運動量方程式の弱形式から、または全ポテンシャルエネルギーの勾配として導出することができる。後者(エネルギー勾配)の方が単純であり、明確に定義されたポテンシャルエネルギーが存在する限り、しばしば用いられる。ここでは超弾性固体に焦点を当てる。

以前述べたように、変形勾配に基づく超弾性エネルギー密度を仮定する。すると、全弾性ポテンシャルエネルギーは次のようになる。

\begin{equation} e = \sum_pV_p^0\Psi_p(\bm{F}_p) \end{equation}

ここで、 V_p^0 は粒子 p の材料空間での体積である。ノード弾性力は、ノード位置で評価された全ポテンシャルエネルギーの負の勾配である。式(186)を用いると、応力ベースの力のMPM空間離散化は次のように与えられる。

\begin{equation} \mathbf{f}_{\bm{i}}(\bm{x}) = -\frac{\partial e}{\partial \bm{x}_{\bm{i}}}(\bm{x}) = -\sum_pV_p^0\left(\frac{\partial \Psi_p}{\partial \bm{F}}(\bm{F}_p(\bm{x}))\right)(\bm{F}_p^n)^T\nabla w_{\bm{i}p}^n \end{equation}

これは、近傍の粒子の弾性応力から生じる、グリッドノード \bm{i} に作用する力である。シンプレクティック・オイラーの場合、力は容易に次のように書ける。

\begin{equation} \mathbf{f}_{\bm{i}}^n = \mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^n) = -\sum_pV_p^0\left(\frac{\partial \Psi_p}{\partial \bm{F}}(\bm{F}_p^n)\right)(\bm{F}_p^n)^T\nabla w_{\bm{i}p}^n \end{equation}

これは、既存の粒子/グリッド重みと粒子の属性に完全に依存する。あるいは、エネルギー密度がない場合、力はコーシー応力を用いて次のように書くこともできる。

\begin{equation} \mathbf{f}_{\bm{i}}^n = \mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^n) = -\sum_pV_p^n\sigma_p^n\nabla w_{\bm{i}p}^n \end{equation}

これは弱形式から来ている。

10-5 MPMスキーム: 完全なアルゴリズム(MPM Scheme: Full Algorithm)

ここでは、シンプレクティック・オイラーMPMの完全な更新スキームの概要を説明する。粒子は既に初期化されていると仮定する(すなわち、各粒子は質量 m_p 、体積 V_p 、初期位置 \bm{x}_p 、初期速度 \bm{v}_p 、アフィン行列 \bm{B}_p 、およびそれらの材料関連パラメータ(例えばネオ・フック材料の場合の \mu\lambda )を持っている)。

  1. 粒子からグリッドへの転送(P2G: Particle-to-Grid Transfer)
    セクション10.1のAPICの式を用いて、最初のステップは粒子の量をグリッドに転送することである。特に、このステップではグリッドの質量と運動量を計算する。
  2. **グリッド速度の計算
    \bm{v}_{\bm{i}} = \frac{m_{\bm{i}}v_{\bm{i}}}{m_{\bm{i}}}。**質量が0のノードについては、 m_{\bm{i}}\bm{v}_{\bm{i}} は手動で0にリセットされる。従来、このステップを実装する際には浮動小数点数の閾値が必要だと考えられていたが、実際には0.0fと直接比較しても不安定な挙動やオーバーフローは生じない。これは、たとえ質量の小さなノードに大きな力がかかっても、その影響は粒子に戻す際に再び重みでスケーリングされるためである。実際、閾値をハードコーディングすると、全運動量のドリフトなど、さまざまな問題を引き起こす。
  3. グリッド自由度の特定
    このステップは実装効率にとって重要である。ゼロでない質量を持つグリッドノードを実際の自由度としてラベル付けを行う。他のすべてのノードは静的なままであり、ソルバーの未知数の一部とは見なされない。
  4. 陽的なグリッド力の計算
    式(189) を用いて、陽的なグリッド力 \mathbf{f}_i^n を計算する。
  5. グリッド速度の更新
    式(183) を用いてグリッド速度を更新する。このステップでは、境界条件や衝突オブジェクトを考慮するべきである。陽的時間積分の場合、各ノード速度はディリクレ境界条件や剛体オブジェクトとの衝突により、独立して所望の値に設定することができる。衝突オブジェクトの扱いに関する詳細はセクション12.1を参照してください。
  6. 粒子変形勾配の更新
    式(181) を用いて粒子の変形勾配を更新する。実際にグリッドを動かしたり、新しいグリッド座標 \bm{x} を計算したりはしないことに注意する。運動は仮想的なものであり、速度のみが陽に保存される。
  7. グリッドから粒子への転送(G2P: Grid-to-Particle Transfer)
    このステップでは、セクション10.1で与えられたスキームを用いて、新しい粒子速度 \bm{v}_p^{n+1} とアフィン行列 \bm{B}_p^{n+1} を計算する。
  8. 粒子の移流(Particle Advection)
    ****最後に、粒子は新しい速度で移流されます: \bm{x}_p^{n+1} = \bm{x}_p^n + \Delta t\bm{v}_p^{n+1} これはAPICが使用される場合にのみ真であることに注意する。FLIPまたはFLIP-PIC混合転送スキームでは、代わりに \bm{x}_p^{n+1} = \Sigma_{\bm{i}}\bm{x}_{\bm{i}}^{n+1}w_{\bm{i}p} を使用すべきである(これら2つの式はPICまたはAPICの場合は同等であることに注意する)。

11 (IMPLICIT TIME INTEGRATION)

陰的時間積分(Implicit time integration)では、陽的時間積分(explicit)との主な違いは、MPMアルゴリズムにおけるグリッド速度の更新ステップである。

11-1 力の微分(Force Derivative)

全ポテンシャルエネルギーは、エネルギー密度 \Psi を用いて次のように表現できることを思い出す。

\begin{equation} \int_{\Omega^0}\Psi(\bm{F}(\bm{X}))\mathrm{d\bm{X}} \end{equation}

ここで \Omega^0 は材料の初期形状である。応力ベースの力のMPM空間離散化は、このエネルギーの離散近似をオイラーグリッドノードの材料位置に関して微分することと等価である。しかし、実際にはオイラーグリッドを変形させないので、グリッドノード位置の変化はグリッドノード速度によって決まると考えることができる。すなわち、もし \bm{x}_{\bm{i}} がグリッドノード \bm{i} の位置であれば、 \hat{\bm{x}}_{\bm{i}} = \bm{x}_{\bm{i}} + \Delta t\bm{v}_{\bm{i}} が、そのノードの現在の速度 \bm{v}_{\bm{i}} が与えられた場合の、そのグリッドノードの変形後の位置となる。すべてのグリッドノード \hat{\bm{x}}_{\bm{i}} のベクトルを \hat{\bm{x}} とすると、全弾性ポテンシャルエネルギーのMPM近似は次のように書ける。

\begin{equation} e(\hat{\bm{x}}) = \sum_pV_p^0\Psi(\hat{\bm{F}}(\hat{\bm{x}})) \end{equation}

ここで V_p^0 は粒子 p が初期に占めていた材料の体積であり、 \hat{\bm{F}}_{Ep} は次のように \hat{\bm{x}} と関連している。

\begin{equation} \hat{\bm{F}}_p(\hat{\bm{x}}) = \left(\bm{I} + \sum_{\bm{i}}(\hat{\bm{x}}_{\bm{i}}-\bm{x}_{\bm{i}})(\nabla w_{\bm{i}p}^n)^T\right)\bm{F}_p^n \end{equation}

この慣例を用いると、応力ベースの力のMPM空間離散化は次のように与えられる。

\begin{equation} -\mathbf{f}_{\bm{i}}(\hat{\bm{x}}) = \frac{\partial e}{\partial \hat{\bm{x}}}_{\bm{i}}(\hat{\bm{x}}) = \sum_pV_p^0\frac{\partial \Psi}{\partial \bm{F}}(\hat{\bm{F}}_p(\hat{\bm{x}}))(\bm{F}_p^n)^T\nabla w_{\bm{i}p}^n \end{equation}

すなわち、 \mathbf{f}_{\bm{i}}(\hat{\bm{x}}) は弾性応力から生じるグリッドノード \bm{i} 上の力である。これはコーシー応力 \bm{\sigma}_p = \frac{1}{J_p^n}\frac{\partial \Psi}{\partial \bm{F}}(\hat{\bm{F}}_p(\hat{\bm{x}}))(\bm{F}_p^n)^T を用いても書ける。[訳注: 繰り返しになりますが、この論文中の \bm{\sigma}_p の定義は一般的ではありません by Gemini]。

\begin{equation} \mathbf{f}_{\bm{i}}(\hat{\bm{x}}) = - \sum_pV_p^n\bm{\sigma}_p\nabla w_{\bm{i}p}^n \end{equation}

ここで、 V_p^n = J_p^nV_p^0 は時刻 t^n で粒子 p が占める材料の体積である。

MPM空間離散化と弾性ポテンシャルのこの関係を強調するのは、グリッド速度 \bm{v}_{\bm{i}} を時間的に陰的に発展させたいからである。この慣例を用いると、 \hat{\bm{x}} に関するポテンシャルのヘッセ行列(Hessian)を利用することにより、更新の弾性部分に対して陰的なステップを取ることができる。このヘッセ行列の任意の増分 \delta\bm{u} に対する作用は、次のように表現できる。

\begin{equation} -\delta\mathbf{f}_{\bm{i}} = \sum_{\bm{j}}\frac{\partial^2e}{\partial\hat{\bm{x}}_{\bm{i}}\partial\hat{\bm{x}}_{\bm{j}}}(\hat{\bm{x}})\delta\bm{u}_{\bm{j}} = \sum_pV_p^0\bm{A}_p(\bm{F}_p^n)^T\nabla w_{\bm{i}p}^n \end{equation}

ここで、

\begin{equation} \bm{A}_p = \frac{\partial^2\Psi}{\partial\bm{F}\partial\bm{F}}(\hat{\bm{F}}_p(\hat{\bm{x}})):\left(\sum_{\bm{j}}\delta\bm{u}_{\bm{j}}(\nabla w_{\bm{j}p}^n)^T\bm{F}_p^n\right) \end{equation}

であり、表記 \bm{A} = \mathcal{C}:\bm{D}
は、 A_{ij} = \mathcal{C}_{ijkl}D_{kl} を意味し、添字 kl については合計が暗黙的に取られる。

これは添字表記でも導出できる。力が式(188)で既に \bm{x} の関数として書かれていることを思い出すと、暗黙的に合計される添字表記で微分すると次のようになる。

\begin{align} \frac{\partial e}{\partial x_{i\alpha}\partial x_{j\tau}} &= -\frac{\partial f_{i\alpha}}{\partial x_{j\tau}} \nonumber\\ &= \sum_pV_p^0\frac{\partial^2\Psi}{\partial F_{\alpha\beta}\partial F_{\tau\alpha}}(\nabla w_{\bm{j}p}^n)_\omega(\nabla w_{\bm{i}p}^n)_{\gamma}(F_p^n)_{\omega\sigma}(F_p^n)_{\gamma\beta} \end{align}

ここで \frac{\partial^2\Psi}{\partial \bm{F}^2}\Psi(\bm{F}) から導出できる。

11-2 後退オイラーシステム(Backward Euler System)

後退オイラー(Backward Euler)時間積分は、 \bm{v}_{\bm{i}}^{n+1} = \bm{v}_{\bm{i}}^n + \Delta t\mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^n)/m_i(式183)を次のように置き換える。

\begin{equation} \bm{v}_{\bm{i}}^{n+1} = \bm{v}_{\bm{i}}^n + \Delta t\mathbf{f}_{\bm{i}}(\bm{x}_{\bm{i}}^{n+1})/m_{\bm{i}} \end{equation}

すなわち、力はグリッドの運動に陰的に(implicitly)依存する。これにより、より大きな時間ステップとより安定した挙動が可能になる。

運動方程式を並べ替えると、グリッドノード速度に関する非線形方程式系が現れる。

\begin{equation} \bm{h}(\bm{v}^{n+1}) = \bm{M}\bm{v}^{n+1} - \Delta t \mathbf{f}(\bm{x}^n + \Delta t \bm{v}^{n+1}) - \bm{M}\bm{v}^n = 0 \end{equation}

11-3 ニュートン法

式(200) を解くことは、伝統的なニュートン・ラフソン法(Newton-Raphson solver)で行うことができる。まず初期推測値 \bm{v}^{(0)}(例えばゼロまたは前の速度 \bm{v}^n)から開始する。そして、解は次のように反復的に改善されていく。

\begin{equation} \bm{v}^{(i+1)}=\bm{v}^{(i)} - \left(\frac{\partial \bm{h}}{\partial \bm{v}}(\bm{v}^{(i)}\right)^{-1}\bm{h}(\bm{v}^{(i)}) \end{equation}

各ステップで、線形システムはMINRESのようなクリロフ部分空間法(Krylov solver)で解かれる。各ニュートン反復において、 \bm{F}_p は速度に応じて更新される必要があり、ニュートン法の解が得られた後の最終的な \bm{F}_p の更新と同様に、各新しい \bm{F}_p の評価のために古い状態 \bm{F}_p^n を保存しておく必要があることに注意する。

多くの場合、ニュートンステップの1回の反復だけで、すでに陽解法よりも桁違いに大きな時間ステップが可能になる。これは、いわゆる半陰的MPM積分(semi-implicit MPM integration)とほぼ同等である。

11-4 線形化された力(Linearized Force)

弾塑性応答は、オイラーグリッドノードの材料位置 \bm{x}_{\bm{i}} = \bm{x}_{\bm{i}}+Δt\bm{v}_{\bm{i}} から定義されると考える。しかし、前に述べたように、このグリッドを実際に変形させることはない。したがって、 \hat{\bm{x}}_{\bm{i}} = \hat{\bm{x}}({\bm{v}} )\bm{v} によって定義されると考えることができる。これを念頭に置いて、次の表記法を使用する: \mathbf{f}_{\bm{i}}^n = \mathbf{f}_{\bm{i}}(\hat{\bm{x}}(0))\mathbf{f}_{\bm{i}}^{n+1} = \mathbf{f}_{\bm{i}}(\hat{\bm{x}}(\bm{v}^{n+1})) および \frac{\partial^2e^n}{\partial \hat{\bm{x}}_{\bm{i}}\partial \hat{\bm{x}}_{\bm{j}}} = -\frac{\partial \mathbf{f}_{\bm{i}}^n}{\partial \hat{\bm{x}}_{\bm{j}}} = -\frac{\partial \mathbf{f}_{\bm{i}}}{\partial \hat{\bm{x}}_{\bm{j}}}(\hat{\bm{x}}(0))

これらの微分を用いて、陰的な更新式を \bm{v}_{\bm{i}}^{n+1} = \bm{v}_{\bm{i}}^n + \Delta tm_{\bm{i}}^{-1}((1-\beta)\mathbf{f}_{\bm{i}}^n + \beta\mathbf{f}_{\bm{i}}^{n+1})\approx\bm{v}_{\bm{i}}^n+\Delta tm_{\bm{i}}^{-1}(\mathbf{f}_{\bm{i}}^n+\beta\Delta t\Sigma_{\bm{j}}\frac{\partial \mathbf{f}_{\bm{i}}^n}{\partial \hat{\bm{x}}_{\bm{j}}}\bm{v}_{\bm{j}}^{n+1}) を使って形成する。これは、 \bm{v}_{\bm{i}}^{n+1} を解くための(質量)対称なシステムにつながる。

\begin{equation} \sum_{\bm{j}}\left(\bm{I}\delta_{\bm{i}\bm{j}} + \beta\Delta t^2m_{\bm{i}}^{-1}\frac{\partial^2e^n}{\partial\hat{\bm{x}}_{\bm{i}}\partial\hat{\bm{x}}_{\bm{j}}}\right)\bm{v}_{\bm{j}}^{n+1} = \bm{v}_{\bm{i}}^* \end{equation}

ここで右辺は、

\begin{equation} \bm{v}_{\bm{i}}^* = \bm{v}_{\bm{i}}^n + \Delta tm_{\bm{i}}^{-1}\mathbf{f}_{\bm{i}}^n \end{equation}

であり、 \beta は陽解法(\beta = 0)、台形則(\beta=\frac{1}{2})、および後退オイラー法(\beta=1)から選択する。このスキームは、ゼロ速度の初期推測値を用いたニュートン法の1ステップを実行することと等価であることが示せる。各時間ステップで1つの線形システムを解く必要がある。半陰的時間積分MPMスキームの概要を図4に示す。

11-5 最適化ベースの時間積分法(Optimization based Integrator)

伝統的なニュートン・ラフソン法(Newton-Raphson solver)は、より標準的な陽解法処理に比べて時間ステップの大幅な改善をもたらしますが、実際には安定性を保つために依然として小さな時間ステップが必要である。

実際、力がポテンシャルエネルギー関数から導かれる場合(ほとんどの場合でそうである)、方程式を解く問題を最適化問題として再定式化することにより、CFL条件程度の任意の \Delta t の下でニュートン法を確実に収束させることが可能です。この最適化問題に対しては、ロバストで効率的な手法が存在する。

次の方程式を解くこと、

\begin{equation} \bm{h}(\bm{v}^{n+1}) = \bm{M}\bm{v}^{n+1} - \Delta t\mathbf{f}(\bm{x}^n+\Delta t\bm{v}^{n+1})-\bm{M}\bm{v}^n = 0 \end{equation}

が以下の目的関数を最小化することと等価であることを示す。

\begin{equation} E(\bm{v}_{\bm{i}}) = \sum_{\bm{i}}\frac{1}{2}m_{\bm{i}}||\bm{v}_{\bm{i}} -\bm{v}_{\bm{i}}^n||^2 + e(\bm{x}_{\bm{i}}^n + \Delta t\bm{v}_{\bm{i}}) \end{equation}

この定式化を用いることで、大幅に大きな時間ステップを取ることが可能になり、元のアプローチへの最小限の変更で大幅な計算節約がもたらされる。この積分器の実装に関する詳細は、[Gast et al., 2015] を参照。

12 その他のトピック(MORE TOPICS)

12-1 衝突オブジェクト(Collision Objects)

グリッド速度に力が適用された直後に、グリッド速度 \bm{v}_{\bm{i}} に対する衝突処理を行う。半陰的時間積分の場合、これは線形システムの右辺に寄与し、衝突しているグリッドノードに対応する自由度は解法の過程で射影除去される。補間による粒子速度とグリッド速度のわずかな不一致を考慮するために、位置を更新する直前に、粒子速度 \bm{v}_p^{n+1} に衝突処理を再度適用することもできる。いずれの場合も、衝突処理は同じ方法で行われる。我々の衝突はすべて非弾性(inelastic)である。粒子ベースの衝突は、粒子の変形勾配に不整合を引き起こすことに注意する。特定のシミュレーションでアーティファクトを生成する可能性があり、必要な場合にのみ有効にするべきである。

衝突オブジェクトはレベルセット(level sets)として表現され、これにより衝突検出(\phi\le0)が簡単になる。衝突の場合、局所的な法線 \bm{n} = \nabla\phi とオブジェクト速度 \bm{v}_{co} が計算される。まず、粒子/グリッド速度 \bm{v} は衝突オブジェクトの参照フレームに変換され、 \bm{v}_{rel} = \bm{v} - \bm{v}_{co} となる。もし物体が離れている場合(v_n = \bm{v}_{rel}\cdot\bm{n}\ge0)、衝突は適用されない。 \bm{v}_t = \bm{v}_{rel}-\bm{n}v_n を相対速度の接線成分とする。固着するインパルスが必要な場合(||\bm{v}_t||\le-\mu v_n)、単純に \bm{v}_{rel}' = 0 とする。ここでプライムは衝突が適用されたことを示す。それ以外の場合は動摩擦を適用し、 \bm{v}_{rel}' = \bm{v}_t + \mu v_n\bm{v}_t/||\bm{v}_t|| となる。ここで \mu は摩擦係数である。最後に、衝突後の相対速度を \bm{v}' = \bm{v}'_{rel} + \bm{v}_{co} で世界座標を戻す。

剛体(rigid)と変形体(deforming)の2種類の衝突オブジェクトを使用した。剛体の場合、静止したレベルセットと、時間とともに変化する可能性のある剛体変換を保存し、これらを用いて任意の点における \phi\bm{n}、および \bm{v}_{co} を計算できる。変形体の場合、レベルセットのキーフレームをロードし、 \phi(\bm{x}, t + \gamma\Delta t) = (1-\gamma)\phi(\bm{x}-\gamma\Delta t\bm{v}_{co}, t)+\gamma\phi(\bm{x}+(1-\gamma)\Delta t\bm{v}_{co}, t+\Delta t) を用いて補間する。ただし、速度は平均速度の代わりに \bm{v}_{co} = (1-\gamma)\bm{v}(\bm{x},t)+\gamma\bm{v}(\bm{x},t+\Delta t) として計算する。

最後に、材料が垂直な面や張り出した面に固着してほしい状況では、一種の固着衝突(sticky collision)を利用する。この場合、クーロン摩擦は不十分である。なぜなら、法線方向の相対速度がゼロ(垂直面)または正(張り出し面で重力により離れようとする)になるためである。

これらの面に対する衝突では、無条件に \bm{v}_{rel}' = 0 とすることでこの効果を実現する。グリッドノードに対するディリクレ境界条件は、固着衝突と等価である。

12-2 ラグランジュ力(Lagrangian Forces)

伝統的なMPMシミュレーションでは、力は式(188) のように計算される。各粒子に変形勾配を保存させることは、容易なトポロジー変化といったメッシュフリー法の利点を持っている。MPMを伝統的なゴムのような弾性固体のシミュレーションに使用する場合、メッシュベースのFEMソルバーのように接続情報を持つ方が理にかなっている。なぜなら、メッシュはより正確な変形勾配計算とより簡単なレンダリングを提供できるからである。ここでは、トポロジー変化を意図しないオブジェクトに対して、メッシュ化されたアプローチも利用可能であり、MPMフレームワークに容易に統合できることを示す。

この場合、全ポテンシャルエネルギー W(\bm{x}_p) を書き下すことができる任意のラグランジュ力モデル(バネ、有限要素など)を使用できる。このラグランジュ力モデルに対応して、力 \mathbf{f}_p=-\frac{\partial W}{\partial\bm{x}_p} を計算する。また、粒子状の任意のベクトル \delta \mu_q が与えられた場合、力の微分 \frac{\partial\mathbf{f}_p}{\partial\bm{x}_q} を掛けて \delta\mathbf{f}_p=\Sigma_q\frac{\partial \mathbf{f}_p}{\partial\bm{x}_q}\delta\bm{\mu}_q を得ることができると仮定する。これらの力関連の構成要素は純粋にラグランジュ的であり、通常の方法で計算されるため、ここでは詳しく説明しない。

メッシュベースの力をラグランジュ力として定義したが、依然としてグリッドを介してそれらを適用する必要がある。力を評価できるように、粒子位置 \bm{x}_p が(概念的に)移動するグリッドノード \bm{x}_{\bm{i}} とどのように関連しているかを記述しなければならない。次に、 \mathbf{f}_p から f_{\bm{i}} を計算し、 \delta\bm{u}_{\bm{i}} が与えられた場合に \delta\mathbf{f}_{\bm{i}} を計算する手段を見つけなければならない。これを行うことで、ラグランジュ力をオイラー力として使用することができる。 \bm{x}_p\bm{x}_{\bm{i}} の更新規則を比較すると、 \bm{x}_p=\Sigma_iw_{\bm{i}p}^n\bm{x}_{\bm{i}} であることがわかる。連鎖律を用いると、

\begin{equation} \mathbf{f}_{\bm{i}} = \sum_pw_{\bm{i}p}^n\mathbf{f}_p\quad\delta\mathbf{f}_{\bm{i}} = \sum_{p,q,\bm{j}}w_{\bm{i}p}^n\frac{\partial\mathbf{f}_p}{\partial\bm{x}_q}w_{\bm{j}q}^n\delta\bm{u}_{\bm{j}} \end{equation}

この式は3つのネストした和で書かれていますが、和を連続して計算することで効率的に計算できる。これらの力はグリッドに適用されるため、MPMアプローチとラグランジュアプローチの両方を同じシミュレーションで使用できる。各粒子は、MPM粒子またはメッシュ化された粒子としてラベル付けされる。メッシュ化された粒子に保存されている変形勾配 \bm{F}_p^n は、それらの粒子についてはこの量がメッシュを使用して計算されるため、決して使用されないことに注意する。これは、MPMをメッシュベースのアプローチと結合するための効果的な手段を提供する。これにより、ラグランジュ法の正確な表面追跡と、オイラーグリッドの自動的な衝突処理が組み合わされる。

Discussion