1. 序論 ― 連続体力学を最適化問題として解く意義
未知ベクトル $\mathbf z$ の最小化問題
\min_{\mathbf z}\; \Phi(\mathbf z)
に対し、ニュートン法は
\Delta\mathbf z
= -\bigl[\nabla^{2}\Phi(\mathbf z_n)\bigr]^{-1}
\nabla\Phi(\mathbf z_n), \qquad
\mathbf z_{n+1}=\mathbf z_n+\Delta\mathbf z
という2次収束アルゴリズムを提供する。本稿の核心は、物理的な平衡方程式と力学的釣り合いを、この形式の連続最適化問題として定式化できることを示すことにある。
この変換視点により、連続最適化の専門家は自身の専門知識を活かして、物理学や工学の難問に直接取り組めるようになる。特に注目すべきは、複雑な物理現象を扱う際に最適化アルゴリズムの洗練された実装をそのまま利用できる点である。
2. 取り扱う変数と仮定
| 記号 | 定義 | 役割 |
|---|---|---|
| $\Omega_0$ | 静止した 3 次元領域 | 初期配置の積分領域 |
| $\mathbf X$ | 物質点の初期位置ベクトル | 参照座標 |
| $\boldsymbol\chi:\Omega_0!\to!\mathbb R^3$ | 変形写像 | 初期位置から変形後位置への写像 |
| $\mathbf u = \boldsymbol\chi - \mathbf X$ | 変位ベクトル場 | 最適化変数 |
| $\mathbf F = \partial\boldsymbol\chi/\partial\mathbf X$ | 変形勾配テンソル | 局所変形を表現する2階テンソル |
| $\mathbf E = \tfrac12(\mathbf F^{!T}\mathbf F - \mathbf I)$ | Green-Lagrange歪テンソル | 有限変形における歪み測度 |
エネルギー密度関数 $\Psi = \Psi(\mathbf F)$ は、変形勾配 $\mathbf F$ の関数として材料の力学的挙動を記述する。
補足 A ― 材料モデルと枠の無関係性
- 材料は履歴を持たない超弾性体とし，エネルギー密度 $\Psi = \Psi(\mathbf F)$ のみで記述する。
- 剛体回転 $\mathbf R!\in!\mathrm{SO}(3)$ に対して $\Psi(\mathbf R\mathbf F) = \Psi(\mathbf F)$ を仮定する（objectivity）。
3. 無限次元を有限次元へ写像する手順
3.1 形状関数による空間離散化
有限要素法では、複雑な連続関数を有限個の節点値と補間関数で近似する。この近似が成立する理論的根拠は次の2点にある：
- Sobolev空間の稠密性: 多項式関数の集合は、十分滑らかな関数空間において稠密である（任意の精度で近似できる）
- 変分原理: エネルギー最小化問題の解は、任意の試験関数に対する仮想仕事の釣り合いと等価である
この理論に基づき、連続的な変位場 $\mathbf u(\mathbf X)$ を離散化するため、形状関数 $N_I(\mathbf X)$ を用いる：
\mathbf u(\mathbf X) = \sum_{I=1}^{n_{\text{dof}}} N_I(\mathbf X)\,\mathbf u_I,
ここで $N_I(\mathbf X)$ は位置 $\mathbf X$ における $I$ 番目の節点の形状関数であり、$\mathbf u_I$ は対応する節点変位ベクトルである。
実際の解析では、以下のような形状関数がよく使われます：
-
1次（線形）要素: 節点でのみ値を持ち、要素内は線形に補間
N_I(\xi) = \frac{1}{2}(1 + \xi_I \xi) \quad \text{（1次元の場合, } \xi \in [-1,1] \text{）} -
2次（2次）要素: 節点と中間点に値を持ち、2次多項式で補間
N_I(\xi) = \begin{cases} \xi(\xi-1)/2 & \text{（左端節点）} \\ (1+\xi)(1-\xi) & \text{（中間節点）} \\ \xi(\xi+1)/2 & \text{（右端節点）} \end{cases} -
Hermite要素: 節点で値と導関数を両方指定し、より滑らかな補間を提供
形状関数は以下の性質を満たす必要があります：
- 節点 $I$ において $N_I(\mathbf X_J) = \delta_{IJ}$ (Kroneckerのデルタ)
- パーティション・オブ・ユニティ: すべての点で $\sum_I N_I(\mathbf X) = 1$
- コンパクトサポート: 各形状関数は局所的にのみ非ゼロ値をとる
これらの性質により、境界条件の適用が容易になり、疎行列計算に適した数値構造が得られます。
これにより無限次元の変位場を有限次元の係数ベクトル
\mathbf u = \{\mathbf u_I\}\in\mathbb R^{n_{\text{dof}}}
へ縮約する。ここで $n_{\text{dof}}$ は自由度の総数である。
形状関数の空間偏微分 $\partial N_I/\partial \mathbf X$ は、有限要素法における変形勾配の計算において重要な役割を果たす。具体的には、変位場 $\mathbf u(\mathbf X)$ の空間微分を次のように計算できる：
\frac{\partial \mathbf u}{\partial \mathbf X} = \sum_{I=1}^{n_{\text{dof}}} \frac{\partial N_I}{\partial \mathbf X} \otimes \mathbf u_I
ここで $\otimes$ はテンソル積を表す。この変位勾配から変形勾配 $\mathbf F$ が次のように計算される：
\mathbf F = \mathbf I + \frac{\partial \mathbf u}{\partial \mathbf X} = \mathbf I + \sum_{I=1}^{n_{\text{dof}}} \frac{\partial N_I}{\partial \mathbf X} \otimes \mathbf u_I
この関係より、変形勾配 $\mathbf F$ と節点変位 $\mathbf u_I$ の関係性が明らかになる。
3.2 ガウス–ルジャンドル積分による数値積分
有限要素法において、要素内の積分を解析的に計算することは一般に困難である。この問題に対処するため、ガウス–ルジャンドル積分が広く採用されている。この方法が有効な理論的根拠は：
- 最適性: $n$ 点のガウス積分は最大 $2n-1$ 次の多項式を正確に積分できる
- 効率性: 最小の評価点数で最高次数の多項式を正確に積分できる
- 誤差の制御: 被積分関数が滑らかであれば、積分点数を増やすことで任意の精度を達成できる
ルジャンドル多項式と積分点の決定
ガウス–ルジャンドル積分の積分点と重みは、ルジャンドル多項式の性質から導出されます。$n$次のルジャンドル多項式 $P_n(x)$ は次の漸化式で定義されます：
P_0(x) = 1, \quad P_1(x) = x
(n+1)P_{n+1}(x) = (2n+1)xP_n(x) - nP_{n-1}(x)
これらの特性により、連続体上の積分を次の近似式で効率的に評価できます：
\int_{\Omega_e} f\,dV
\;\approx\;
\sum_{g=1}^{n_g} w_g\,f(\boldsymbol\xi_g).
ここで $\Omega_e$ は要素領域、$w_g$ は重み係数、$\boldsymbol\xi_g$ は積分点（ガウス点）の位置、$n_g$ は積分点数である。
実際の計算では、まず要素を標準領域（例：[-1,1]^3の立方体）に写像し、その上でガウス積分を適用します。一般的な積分点と重みの組み合わせは：
1次元の例：
- 1点積分: $\xi_1 = 0$, $w_1 = 2$ (1次多項式まで正確)
- 2点積分: $\xi_{1,2} = \mp\frac{1}{\sqrt{3}}$, $w_{1,2} = 1$ (3次多項式まで正確)
- 3点積分: $\xi_1 = 0$, $\xi_{2,3} = \pm\sqrt{\frac{3}{5}}$, $w_1 = \frac{8}{9}$, $w_{2,3} = \frac{5}{9}$ (5次多項式まで正確)
多次元の場合は、1次元ガウス積分のテンソル積として定義されます：
\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1} f(\xi,\eta,\zeta)\,d\xi\,d\eta\,d\zeta \approx
\sum_{i=1}^{n_i}\sum_{j=1}^{n_j}\sum_{k=1}^{n_k} w_i w_j w_k\,f(\xi_i,\eta_j,\zeta_k)
積分点の選択は精度と計算効率のトレードオフに直接影響します：
- 積分点が少なすぎると、特に高次要素や非線形問題で不正確な結果を生じる
- 積分点が多すぎると、計算コストが増加する
- 特定の問題（例：非圧縮性材料）では、選択的低減積分が数値安定性を向上させる
この数値積分法により、複雑な形状や非線形材料特性を持つ問題でも効率的に高精度な解が得られます。
補足 B ― 積分点数とロッキング
次数 $p$ の形状関数には
$n_g \ge \lceil (2p+1)/2\rceil^3$
が経験的最小点数である。一次要素で大変形解析を行う場合、体積ロッキングを抑制する目的で選択積分や混合変数法の導入が推奨される。
4. 全ポテンシャルエネルギーとニュートン法
4.1 連続体としてのエネルギー
連続体の全ポテンシャルエネルギー $\Pi[\mathbf u]$ は変分原理に基づく。この原理は、平衡状態にある系では全ポテンシャルエネルギーが停留値（最小値）をとるという物理法則に由来する。具体的には：
- エネルギー保存則: 閉じた系のエネルギーは保存される
- 最小作用の原理: 自然界の過程は作用積分を最小にする経路で進む
- 仮想仕事の原理: 平衡状態では仮想変位に対する仮想仕事の総和がゼロになる
これらの物理原理から、連続体の全ポテンシャルエネルギーは次式で表される：
\Pi[\mathbf u] =
\int_{\Omega_0} \Psi(\mathbf F)\,dV
- \int_{\Omega_0} \mathbf f_B\!\cdot\!\mathbf u\,dV
- \int_{\Gamma_0^t} \mathbf f_S\!\cdot\!\mathbf u\,dA.
この式の各項は以下の物理的意味を持ちます：
-
第1項: 内部エネルギー（歪みエネルギー）
- 材料内部に蓄えられるエネルギーで、変形に抵抗する力の源
- $\Psi(\mathbf F)$ は単位体積あたりの歪みエネルギー密度関数
-
第2項: 体積力によるポテンシャルエネルギー
- 物体全体に作用する力（例：重力）による仕事の負値
- $\mathbf f_B$ は単位体積あたりの体積力ベクトル
-
第3項: 表面力によるポテンシャルエネルギー
- 物体表面に作用する力（例：圧力）による仕事の負値
- $\mathbf f_S$ は単位面積あたりの表面力ベクトル
- $\Gamma_0^t$ は力が作用する境界面
全ポテンシャルエネルギー最小化の原理により、平衡状態では $\Pi[\mathbf u]$ の第一変分がゼロになります：
\delta\Pi[\mathbf u] = 0 \quad \forall \delta\mathbf u
これを展開すると：
\int_{\Omega_0} \frac{\partial\Psi}{\partial\mathbf F}:\delta\mathbf F\,dV
- \int_{\Omega_0} \mathbf f_B\!\cdot\!\delta\mathbf u\,dV
- \int_{\Gamma_0^t} \mathbf f_S\!\cdot\!\delta\mathbf u\,dA = 0
この式は、力の釣り合い（運動量保存則）のポテンシャル形式です。したがって、全ポテンシャルエネルギーの停留値を求めることは、連続体の力学的平衡状態を求めることと等価になります。
この変分原理が有限要素法の基礎となり、次節で示す離散化されたエネルギー関数の最小化問題へと帰着します。
4.2 離散スカラー関数
空間離散化と数値積分により、連続体のエネルギーは有限次元のスカラー関数に変換される：
\Phi(\mathbf u)=
\sum_e \sum_g w_g\,
\Psi\bigl(\mathbf F(\boldsymbol\xi_g;\mathbf u)\bigr)
-
\mathbf u^{\mathsf T}\mathbf f_{\text{ext}}.
ここで $\mathbf f_{\text{ext}}$ は外力の離散表現であり、先の体積力と表面力を反映したものである。
4.3 ニュートン法による最適化 ― 連続最適化理論の直接適用
離散化されたエネルギー $\Phi(\mathbf u)$ の最小化は、連続最適化問題そのものである。ここで、最適化理論のあらゆる知見が直接適用可能となる。ニュートン法による更新式は：
\Delta\mathbf u
= -\bigl[\nabla^{2}\Phi(\mathbf u_n)\bigr]^{-1}
\nabla\Phi(\mathbf u_n), \qquad
\mathbf u_{n+1}=\mathbf u_n+\Delta\mathbf u.
この更新式は、連続最適化の標準的なニュートン法と完全に一致する。その物理的解釈は：
- $\nabla\Phi(\mathbf u_n)$ は力の不均衡（残差力）
- $\nabla^{2}\Phi(\mathbf u_n)$ は系の剛性（局所的な曲率）
- $\Delta\mathbf u$ は平衡状態への修正量
連続最適化に精通した読者にとって、この対応関係は極めて有益である。以下のような最適化技術がそのまま適用できる：
- 準ニュートン法: BFGS/L-BFGS法でヘッセン行列の計算・保存コストを削減
- 信頼領域法: 大きな変形増分を安定して解く際に有効
- 制約付き最適化: 接触条件や拘束条件を直接組み込める
- 大規模問題の解法: 前処理付き共役勾配法や多重格子法などの手法がそのまま使用可能
- 並列計算戦略: ドメイン分解法やマルチスタート法などの最適化並列化手法
また、現代的な最適化ソルバー（IPOPT, KNITRO, NLoptなど）や自動微分ツールを活用することで、複雑な非線形材料モデルでも効率的に解析できます。
この「最適化としての物理」という視点は、両分野の手法を融合させ、複雑な非線形問題に対する新たなアプローチを可能にします。
補足 C ― 収束対策と分岐検出
- ヘッセ行列が半正定値の場合や荷重増分が大きい場合，反復が発散し得る。線形検索・ダンピング・アークレングス法が有効である。
- $\det\nabla^{2}\Phi \to 0$ は分岐点接近の指標となるため，最小固有値を監視してクリティカルロードを推定できる。
5. 勾配とヘッセ行列の具体形
5.1 残差ベクトル
エネルギー関数 $\Phi(\mathbf u)$ の勾配（残差ベクトル）は：
\mathbf r(\mathbf u)=
\nabla\Phi(\mathbf u) =
\sum_e \sum_g w_g\,
\mathbf B_g^{\mathsf T}\,\mathbf P_g
- \mathbf f_{\text{ext}},
ここで：
- $\mathbf B_g$ は形状関数の微分から構成される行列で、変形勾配 $\mathbf F$ と節点変位 $\mathbf u$ の関係を表現する。具体的には：
\mathbf B_g = \frac{\partial \mathbf F}{\partial \mathbf u} =
\begin{bmatrix}
\frac{\partial \mathbf F}{\partial \mathbf u_1} & \frac{\partial \mathbf F}{\partial \mathbf u_2} & \cdots & \frac{\partial \mathbf F}{\partial \mathbf u_{n_{\text{dof}}}}
\end{bmatrix}
各成分は形状関数の空間微分から直接導出される：
\frac{\partial \mathbf F}{\partial \mathbf u_I} = \frac{\partial}{\partial \mathbf u_I}\left(\mathbf I + \sum_{J=1}^{n_{\text{dof}}} \frac{\partial N_J}{\partial \mathbf X} \otimes \mathbf u_J\right) = \frac{\partial N_I}{\partial \mathbf X} \otimes \mathbf I
ここで $\mathbf I$ は単位テンソルである。このように、$\mathbf B_g$ 行列の各成分は形状関数 $N_I$ の空間微分から直接構成される。
- $\mathbf P_g$ は積分点での第1 Piola-Kirchhoff応力テンソルで、次式で定義される：
\mathbf P = \frac{\partial\Psi}{\partial\mathbf F}\bigl(\mathbf F\bigr)
したがって、残差ベクトルは基本的に各節点での変位に対するエネルギーの変化率を表し、形状関数の空間微分を通じて計算される。
5.2 一貫接線剛性
エネルギー関数 $\Phi(\mathbf u)$ のヘッセ行列（接線剛性行列）は：
\mathbf K(\mathbf u) =
\nabla^{2}\Phi(\mathbf u) =
\sum_e \sum_g w_g\,
\mathbf B_g^{\mathsf T}\,
\mathbf D_g\,
\mathbf B_g\,
|\det\mathbf J_g|,
ここで：
- $\mathbf D_g$ は材料剛性テンソルで、次式で定義される：
\mathbf D = \frac{\partial^2\Psi}{\partial\mathbf F^2}\bigl(\mathbf F\bigr)
- $\mathbf J_g$ はヤコビ行列で、参照座標から実座標への変換を表す
$\mathbf K$ はニュートン反復におけるヘッセ行列に相当し、完全微分を用いることで 2 次収束が保証される。
6. 小変形理論への極限
変位が微小なとき、次の近似が成り立つ：
\mathbf F \approx \mathbf I + \nabla\mathbf u, \quad
\Psi(\mathbf F) \approx \tfrac12\,\mathbf E:\mathbf C:\mathbf E,
ここで $\mathbf C$ は材料の弾性係数テンソルである。これを用いると、要素剛性は：
\mathbf K_e =
\int_{\Omega_e}
\mathbf B_L^{\mathsf T}\,
\mathbf C\,
\mathbf B_L
\,d\Omega.
ここで $\mathbf B_L$ は線形歪み-変位行列であり、変位の勾配と線形歪みテンソルを関連付ける。これは線形有限要素法における剛性行列の標準形式に一致する。
6.1 Voigt記法と教科書的表現
小変形理論では、テンソル表記からマトリックス表記への変換がよく行われる。具体的にはVoigt記法を用いて：
\mathbf C_{ijkl} \rightarrow \mathbb{C}_{IJ}, \quad
\mathbf E_{ij} \rightarrow \boldsymbol{\varepsilon}_I
ここで添字変換は $(11) \rightarrow 1, (22) \rightarrow 2, (33) \rightarrow 3, (23) \rightarrow 4, (13) \rightarrow 5, (12) \rightarrow 6$ である。
この記法を用いると、線形弾性問題の標準的な表現として次の関係式が得られる：
\boldsymbol{\sigma} = \mathbb{C} \boldsymbol{\varepsilon}
ここで $\boldsymbol{\sigma}$ と $\boldsymbol{\varepsilon}$ はそれぞれ応力と歪みのVoigt表記ベクトルであり：
\boldsymbol{\sigma} =
\begin{bmatrix}
\sigma_{11} \\
\sigma_{22} \\
\sigma_{33} \\
\sigma_{23} \\
\sigma_{13} \\
\sigma_{12}
\end{bmatrix}, \quad
\boldsymbol{\varepsilon} =
\begin{bmatrix}
\varepsilon_{11} \\
\varepsilon_{22} \\
\varepsilon_{33} \\
2\varepsilon_{23} \\
2\varepsilon_{13} \\
2\varepsilon_{12}
\end{bmatrix}
等方弾性体の場合、剛性マトリックス $\mathbb{C}$ は次のように表される：
\mathbb{C} =
\begin{bmatrix}
\lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\
\lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\
\lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\
0 & 0 & 0 & \mu & 0 & 0 \\
0 & 0 & 0 & 0 & \mu & 0 \\
0 & 0 & 0 & 0 & 0 & \mu
\end{bmatrix}
ここで $\lambda$ と $\mu$ はLaméの弾性定数である。線形歪み-変位マトリックス $\mathbf B_L$ は形状関数の偏微分を用いて構成される：
\mathbf B_L^{(I)} =
\begin{bmatrix}
\frac{\partial N_I}{\partial x} & 0 & 0 \\
0 & \frac{\partial N_I}{\partial y} & 0 \\
0 & 0 & \frac{\partial N_I}{\partial z} \\
0 & \frac{\partial N_I}{\partial z} & \frac{\partial N_I}{\partial y} \\
\frac{\partial N_I}{\partial z} & 0 & \frac{\partial N_I}{\partial x} \\
\frac{\partial N_I}{\partial y} & \frac{\partial N_I}{\partial x} & 0
\end{bmatrix}
これを用いて、要素剛性行列は以下の教科書的な形で表現される：
\mathbf K_e = \int_{\Omega_e} \mathbf B_L^{\mathsf T} \mathbb{C} \mathbf B_L \, d\Omega
また、節点力ベクトルは：
\mathbf f_e = \int_{\Omega_e} \mathbf B_L^{\mathsf T} \boldsymbol{\sigma} \, d\Omega
これらは、構造解析の教科書で通常見られる標準形式である。
Voigt記法の長所と短所
長所:
- 計算効率の向上: テンソル演算を効率的な行列演算に変換
- メモリ使用量の削減: 対称性を活用して冗長な情報を排除
- プログラミングの容易さ: 標準的な線形代数ライブラリで実装可能
短所:
- 物理的直感の低下: テンソルの幾何学的意味が見えにくくなる
- せん断成分の非一貫性: 変換規則が歪みと応力で異なる場合がある
- 一般化の難しさ: 高次テンソルへの拡張が複雑になる
このように、小変形理論ではVoigt記法により、テンソル表記（5節）からマトリックス表記（6節）への変換が行われる。この変換は実装上の効率性をもたらす一方で、非線形理論と線形理論の間で記号体系の一貫性を損なう要因となる。
7. まとめ
- 形状関数と数値積分により、連続体力学の平衡を有限次元スカラー関数 $\Phi(\mathbf u)$ に帰着した。
- ニュートン法をそのまま適用することで効率的に平衡解を得られる。
- 勾配・ヘッセ行列は、それぞれ応力 $\mathbf P$ と材料剛性 $\mathbf D$ の形で物理的意味をもつ。
- 小変形極限では標準的な線形 FEM への収束が確認できる。