幾何学非線形問題の前に、まずよく用いる微小変形理論を復習しましょう。図 1 は典型的な静力学的な境界値問題を示しています。連続体は領域 $\Omega$ を占め、その境界が $\partial \Omega$ です。境界は変位境界 $\partial \Omega_{\mathbf{u}}$ と力境界 $\partial \Omega_{\mathrm{a}}$ に分割され（ $\partial \Omega_{\mathbf{u}} \cup \partial \Omega_{\mathrm{a}}=\partial \Omega 、 \partial \Omega_{\mathbf{u}} \cap \partial \Omega_{\mathrm{a}}=\varnothing$ ）、変位境界では変位 $\overline{\boldsymbol{u}}$ が与えられ、力境界 $\partial \Omega_{\mathrm{a}}$ では表面力 $\mathbf{f}_{a}$ が与えられます。連続体領域 $\Omega$ には体積力が作用することもあり、 $\mathbf{f}_{v}$ と書きます。未知量は変位場 $\boldsymbol{u}(\boldsymbol{x})$ であり（ $\boldsymbol{x}$ が位置ベクトルとします）、各材料点が変形後にどれだけ移動したかを表します。

$$
\begin{gathered}
& \partial \Omega_{\mathbf{u}} \cup \partial \Omega_{\mathrm{a}}=\partial \Omega \\
& \partial \Omega_{\mathbf{u}} \cap \partial \Omega_{\mathrm{a}}=\varnothing\\ \\
& 図 1．微小変形理論における静力学的な境界値問題
\end{gathered}
$$

　外力（体積力 $\mathbf{f}_{v}$ または表面力 $\mathbf{f}_{a}$ ）が作用すると、構造物が変形します。変形がどれほど複雑でも、各材料点では力のつり合いが満たされなければなりません。慣性項を無視した静力学における線運動量保存（Cauchy の第 1 運動法則）を満たす必要があります。

$$
\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}+\mathbf{f}_{v}=\mathbf{0} \quad\text { in } \quad \Omega $$ $\boldsymbol{\sigma}$ が応力テンソル（2 階テンソル）、 $\boldsymbol{\nabla} \cdot(\cdot)$ は発散を表します。 $\mathbf{f}_{v}$ は単位体積当たりの体積力です。基底ベクトル $\boldsymbol{e}_{j}$ を用いて展開すると、以下のようになります。
$$
\frac{\partial \sigma_{i j}}{\partial x_{i}} \boldsymbol{e}_{j}+\left(\mathrm{f}_{v}\right)_{j} \boldsymbol{e}_{j}=0_{j} \boldsymbol{e}_{j} \quad\text { in } \quad\Omega
$$

また、角運動量のつり合いより $\boldsymbol{\sigma}$ の対称性が得られます（Cauchy の第 2 運動法則： $\boldsymbol{\sigma}=\boldsymbol{\sigma}^{T}$ ）。本記事では、幾何学非線形の影響に焦点を当てるため、構成則に関しては線形弾性体を仮定します。応力とひずみは以下のような線形関係を持つことにします。$\mathbb{C}$ が 4 階テンソルであり、 $\boldsymbol{\varepsilon}$ が微小ひずみテンソルとなります。
$$
\boldsymbol{\sigma}=\mathbb{C}: \boldsymbol{\varepsilon} \quad \text { in } \quad \Omega
$$

　図 2 に示すように、変位の絶対量よりも局所的な変形を表すには、空間における変位の変化率、すなわち変位勾配 $\boldsymbol{\nabla} \boldsymbol{u}$ が重要となります（ $\boldsymbol{\nabla} \boldsymbol{u} \equiv \boldsymbol{u} \otimes \boldsymbol{\nabla}$ ）。図 2 では物体内部点 P の無限小近傍における $\boldsymbol{u}$ の勾配によって生じる 2 つのベクトル $d \boldsymbol{x}_{1}, d \boldsymbol{x}_{2}$ の内積の変化の様子を表しています。つまり、局所的な変形は $\boldsymbol{u}$ そのものではなく $\partial \mathbf{u} / \partial x_{i}$ によって決まります。このため、微小変形理論ではひずみが $\boldsymbol{\nabla} \boldsymbol{u}$ から構成されます。

　変位勾配テンソル $\boldsymbol{\nabla} \boldsymbol{u}$ の対称部分を取ると、以下の対称な微小ひずみテンソル $\boldsymbol{\varepsilon}$（2 階テンソル） が得られます。これより、ひずみ $\boldsymbol{\varepsilon}$ と $\boldsymbol{u}$ の線形関係も確認できます。

$$
\boldsymbol{\varepsilon}=\frac{1}{2}\left[\boldsymbol{\nabla} \boldsymbol{u}+(\boldsymbol{\nabla} \boldsymbol{u})^{T}\right] \quad \text { in } \quad \Omega
$$

また、デカルト座標系を使用すると、このひずみテンソルの成分表示ができるようになります。

$$
\begin{aligned}
\boldsymbol{\varepsilon} & =\frac{1}{2}\left(\boldsymbol{\nabla} \boldsymbol{u}+(\boldsymbol{\nabla} \boldsymbol{u})^{T}\right) \\
& =\frac{1}{2}\left(\frac{\partial \boldsymbol{u}}{\partial x_{j}} \otimes \boldsymbol{e}_{j}+\boldsymbol{e}_{i} \otimes \frac{\partial \boldsymbol{u}}{\partial x_{i}}\right) \\
& =\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right) \boldsymbol{e}_{i} \otimes \boldsymbol{e}_{j}
\end{aligned}
$$

より直観的に見るために、行列の形に書き換えます。
$$
[\boldsymbol{\varepsilon}]=\left[\begin{array}{ccc}
\frac{\partial u_{1}}{\partial x_{1}} & \frac{1}{2}\left(\frac{\partial u_{1}}{\partial x_{2}}+\frac{\partial u_{2}}{\partial x_{1}}\right) & \frac{1}{2}\left(\frac{\partial u_{1}}{\partial x_{3}}+\frac{\partial u_{3}}{\partial x_{1}}\right) \\
\frac{1}{2}\left(\frac{\partial u_{2}}{\partial x_{1}}+\frac{\partial u_{1}}{\partial x_{2}}\right) & \frac{\partial u_{2}}{\partial x_{2}} & \frac{1}{2}\left(\frac{\partial u_{2}}{\partial x_{3}}+\frac{\partial u_{3}}{\partial x_{2}}\right) \\
\frac{1}{2}\left(\frac{\partial u_{3}}{\partial x_{1}}+\frac{\partial u_{1}}{\partial x_{3}}\right) & \frac{1}{2}\left(\frac{\partial u_{3}}{\partial x_{2}}+\frac{\partial u_{2}}{\partial x_{3}}\right) & \frac{\partial u_{3}}{\partial x_{3}}
\end{array}\right]
$$

幾何学的な意味を確認するには、 2 次元問題に落とすことも一つの方法だと考えられます（元の 3 次元問題の理解も損ないません）。
$$
[\boldsymbol{\varepsilon}]=\left[\begin{array}{cc}
\frac{\partial u_{1}}{\partial x_{1}} & \frac{1}{2}\left(\frac{\partial u_{1}}{\partial x_{2}}+\frac{\partial u_{2}}{\partial x_{1}}\right) \\
\frac{1}{2}\left(\frac{\partial u_{2}}{\partial x_{1}}+\frac{\partial u_{1}}{\partial x_{2}}\right) & \frac{\partial u_{2}}{\partial x_{2}}
\end{array}\right]
$$

図 3 は 2 次問題における代表的な変形モードを示し、微小ひずみテンソルの各成分の幾何学的な意味を理解するための例になります。

最初に紹介した2種類の境界条件（力と変位）を併せると、次の境界値問題が得られます。

$$
\begin{array}{r}
\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}+\mathbf{f}_{v}=\mathbf{0} \quad \text { in } \quad \Omega \\ \boldsymbol{\sigma}=\mathbb{C}: \boldsymbol{\varepsilon} \quad \text { in } \quad \Omega \\ \boldsymbol{\varepsilon}=\frac{1}{2}\left[\boldsymbol{\nabla} \boldsymbol{u}+(\boldsymbol{\nabla} \boldsymbol{u})^{T}\right] \quad \text { in } \quad \Omega \\ \boldsymbol{\sigma} \cdot \boldsymbol{n}=\mathbf{f}_{a} \quad \text { on } \quad \partial \Omega_{\mathrm{a}} \\
\boldsymbol{u}=\overline{\boldsymbol{u}} \quad \text { on } \quad \partial \Omega_{\mathbf{u}}
\end{array}
$$

　これがよく使用される微小変形理論の支配方程式であり、強形式とも言われています。この問題は外力から始まり、変形が生じ、各材料点で力と力のモーメントのつり合いが満たされています。しかし、よく考えると素朴な疑問が生まれます。このつり合いは、いったいどの形状（幾何学）に対して成り立っているのでしょうか？つまり、微小変形理論では（これまでの議論）、力のつり合いを考える際に幾何形状の変化の影響を無視していたことが問題でした。
　例えば、幾何形状の変化を含むように模式図を描き直してみると、図 4 が示すようになります。

$$
\begin{gathered}
& \partial \Omega_{\mathrm{u}}^{\mathrm{t}} \cup \partial \Omega_{\mathrm{a}}^{\mathrm{t}}=\partial \Omega \\
& \partial \Omega_{\mathrm{u}}^{\mathrm{t}} \cap \partial \Omega_{\mathrm{a}}^{\mathrm{t}}=\varnothing\\ \\
& 図 4．有限変形理論における静力学的な境界値問題
\end{gathered}
$$

図の灰色領域は、時刻 t における現在配置（変形後の領域）$\Omega_{\mathrm{t}}$ を表します。破線で囲む領域 $\Omega_{0}$ は変形前の形状を表します。幾何学形状の変化を含むように描くと、外力は変更後の形状に作用しないといけないことが明らかになりました。つまり、 $\boldsymbol{n}$ が現在配置 $\Omega_{\mathrm{t}}$ の境界上の外向き単位法線ベクトルであり、 $\mathbf{f}_{v}$ が現在体積で定義される単位体積当たりの力となります。 2 種類の境界条件（変位境界、力境界） も現在配置で定義されるようになります（ $\partial \Omega_{\mathbf{u}}^{t} \cup \partial \Omega_{\mathrm{a}}^{t}=\partial \Omega 、 \partial \Omega_{\mathrm{u}}^{t} \cap \partial \Omega_{\mathrm{a}}^{t}=\varnothing$ ）。

　力のつり合いについて、微小変形（幾何学的線形）理論では、ほぼ不変とみなせる計算領域 $\Omega$上で通常次のように書きますが、
$$
\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}+\mathbf{f}_{v}=\mathbf{0} \quad \text { in } \quad \Omega
$$

大変位／大回転を許すと、変形前後の形状の違いが無視できなくなり、真の力のつり合いは現在配置 $\Omega_{t}$ 上で成立し、発散演算子も $\boldsymbol{x}$ に関して取る必要があります。したがって上式は次のように書き換えられます：
$$
\boldsymbol{\nabla}_{\boldsymbol{x}} \cdot \boldsymbol{\sigma}+\mathbf{f}_{v}=\mathbf{0} \quad \text { in } \quad \Omega_{\mathrm{t}}
$$

ここで $\boldsymbol{\sigma}$ は Cauchy 応カテンソル（真応カテンソル）となります。 $\Omega_{\mathrm{t}}$ は未知の変位場 $\boldsymbol{u}$ によって決まるため、平衡方程式の定義域そのものが未知になります。これは幾何学的非線形が求解難度を大きく上げる理由の一つです。
　ひずみも変形前と変形後を考慮するように定義する必要があります。より見やすくするために、ここでは直接位置ベクトルが存在することを前提として式展開を行います。運動 $\varphi(\mathbf{X}, t)$ を現在配置の位置べクトルにします ：
$$
\begin{gathered}
\boldsymbol{x}=\varphi(\mathbf{X}, t), \quad \boldsymbol{X} \in \Omega_{0}, \quad \boldsymbol{x} \in \Omega_{\mathrm{t}} \\
\boldsymbol{x}(\boldsymbol{X}, t)=\boldsymbol{X}+\boldsymbol{u}(\boldsymbol{X}, t)
\end{gathered}
$$
$\boldsymbol{u}(\boldsymbol{X}, t)$ は参照座標 $\boldsymbol{X}$ を独立変数とする変位場（ラグランジュ記述）とします。 $\boldsymbol{u}$ が大きくても、貫通や自己交差が生じない限り、写像 $\varphi(\mathbf{X}, t)$ が 1 対 1 で保たれ、この表現は有効となります。微小変形理論のときと同じように、局所の変形を測るには、運動 $\varphi(\mathbf{X}, t)$ に注目するよりも、各点での微分つまり $d \boldsymbol{X}$ と $d \boldsymbol{x}$ の関係に注目する方が本質的です。そこで、以下の線形変換を定義します。
$$
\begin{gathered}
d \boldsymbol{x}=\boldsymbol{F} d \boldsymbol{X} \\
d \boldsymbol{X}=\boldsymbol{F}^{-1} d \boldsymbol{x}
\end{gathered}
$$

$\boldsymbol{F}$ は変形勾配テンソルと呼ばれます。一般的には参照配置と現在配置の双方で同じデカルト直交基 $\mathbf{e}_{i}$ を選ぶので、変形勾配テンソルは次のように表すことができます： $$ \boldsymbol{F}=\frac{\partial x^{i}}{\partial X^{J}} \mathbf{e}_{i} \otimes \mathbf{e}_{J} \equiv F_{i J} \mathbf{e}_{i} \otimes \mathbf{e}_{J}
$$

変形勾配テンソル $\boldsymbol{F}$ を用いて任意の 2 つのベクトル $d \boldsymbol{x}_{1}, d \boldsymbol{x}_{2}$ の内積の変化を表してみると、
$$
\begin{aligned}
& d \boldsymbol{x}_{1} \cdot d \boldsymbol{x}_{2}-d \boldsymbol{X}_{1} \cdot d \boldsymbol{X}_{2} \\
& =d \boldsymbol{x}_{1} \cdot d \boldsymbol{x}_{2}-d \boldsymbol{x}_{1} \boldsymbol{F}^{-T} \cdot \boldsymbol{F}^{-1} d \boldsymbol{x}_{2} \\
& =d \boldsymbol{x}_{1} \cdot d \boldsymbol{x}_{2}-d \boldsymbol{x}_{1}\left(\boldsymbol{F F}^{T}\right)^{-1} d \boldsymbol{x}_{2} \\
& =d \boldsymbol{x}_{1} \cdot d \boldsymbol{x}_{2}-d \boldsymbol{x}_{1}\left(\boldsymbol{V}^{2}\right)^{-1} d \boldsymbol{x}_{2} \\
& =d \boldsymbol{x}_{1}\left(\mathbf{1}-\boldsymbol{V}^{-2}\right) d \boldsymbol{x}_{2}
\end{aligned}
$$

ここで、 $\boldsymbol{V}^{2} \equiv \boldsymbol{F} \boldsymbol{F}^{T}$ の関係を利用しました。 $\boldsymbol{V}$ は左ストレッチテンソルと呼ばれています（より詳細な内容は非線形固体力学または連続体力学の専門書をご参照ください）。上式から（ $\mathbf{1}-\boldsymbol{V}^{-2}$ ）がひずみテンソルの役割を持つことがわかるかと思います。つまり、 $\boldsymbol{V}$ を用いると以下の対数型のひずみテンソルを定義することができます。
$$
\boldsymbol{\varepsilon}=\ln \boldsymbol{V}
$$

　このように定義されたひずみは対数ひずみテンソルまたはHencky ひずみテンソルと呼ばれます。 $\boldsymbol{F}$ が変位 $\boldsymbol{u}$ の関数であるため、 $\boldsymbol{V}$ が変位 $\boldsymbol{u}$ の非線形関数であることもわかります。微小変形理論のときと同じように、材料線形を仮定すると、以下の構成関係が得られます。
$$
\boldsymbol{\sigma}=J^{-{1}} \mathbb{C}: \boldsymbol{\varepsilon}, \quad J=\operatorname{det} \boldsymbol{F}
$$

注意していただきたいのは、幾何学非線形性を考慮するとき、材料線形性を仮定しても、Hencky ひずみテンソルと Cauchy 応カテンソルの間には $J^{-1}$ の項が増えています。境界条件を追加して、すべての量を現在配置 $\Omega_{\mathrm{t}}$ 上で記述すると、支配方程式は次のようにまとめられます：
$$
\begin{array}{r}
\boldsymbol{\nabla}_{\boldsymbol{x}} \cdot \boldsymbol{\sigma}+\mathbf{f}_{v}=\mathbf{0} \quad \text { in } \quad \Omega_{\mathrm{t}} \\
\boldsymbol{\sigma}=J^{-{1}} \mathbb{C}: \boldsymbol{\varepsilon} \quad \text { in } \quad \Omega_{\mathrm{t}} \\
\boldsymbol{\varepsilon}=\ln \boldsymbol{V} \quad \text { in } \quad \Omega_{\mathrm{t}} \\
\boldsymbol{\sigma} \cdot \boldsymbol{n}=\mathbf{f}_{\boldsymbol{a}} \quad \text { on } \quad \partial \Omega_{\mathrm{a}}^{\mathrm{t}} \\
\boldsymbol{u}=\overline{\boldsymbol{u}} \quad \text { on } \quad \partial \Omega_{\mathrm{u}}^{\mathrm{t}}
\end{array}
$$

　ここでは材料非線形など他の非線形性を一旦取り除いているため、この支配方程式の非線形性の主な起源は次の 2 点に整理できます：（1）変位とひずみの関係が非線形であること。（2）解くべき領域 $\Omega_{\mathrm{t}}$ 自体が未知であること。このような点に対応するため、COMSOL Multiphysics ${ }^{\text {® }}$ における幾何学的非線形性を考慮した大変形固体力学では Total Lagrange などの形式を用い、弱形式を含む平衡方程式を参照構形へ引き戻して扱います。次の記事で続けて紹介します。

著者

張広志 Ph．D．
計測エンジニアリングシステム株式会社 技術部
非線形性固体力学•非線形有限要素法ならびに、COMSOL Multiphysics ${ }^{\text {® }}$ への技術サポートに従事
計算力学技術者 固体力学分野 1 級 振動分野 1 級