概要
土木の分野において、基本的には微小変形理論の範囲で解析が行われることが多いです。しかし、地盤などの大変形を扱う場合、有限変形理論の枠組みで扱う必要があります。
今回の記事では、微小変形理論で使われる線形弾性モデルを有限変形理論に拡張した超弾性モデル、hencky モデルについて簡単に紹介したいと思います。
hencky ひずみと応力
Piola-Kirchhoff 応力テンソルは、エネルギー関数 $\psi$ を使い
\begin{align}
P = \frac{\partial \psi}{\partial F}
\end{align}
とかける。 hencky ひずみにおけるエネルギー関数 $\psi$ は
\begin{align}
\psi(\epsilon) &= \frac{\lambda}{2} (\mathrm{Tr}\epsilon)^2 + \mu \epsilon\colon\epsilon \\
&=\frac{1}{2} \epsilon\colon D \colon \epsilon
\end{align}
と書ける。$F$ は変形勾配であり、$\epsilon$ はひずみである。また $\colon$ は、具体的に成分で書くと
\begin{align}
\epsilon\colon\epsilon = \epsilon_{ij} \epsilon_{ij}
\end{align}
を意味する。$D$ は弾性テンソルであり、
\begin{align}
D_{ijkl}=G(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}) + \left(K-\frac{2}{3}G \right)\delta_{ij}\delta_{kl}
\end{align}
である。$K$ は体積弾性係数、$G$ はせん断弾性係数である。また、$\lambda,\mu$ はラメ定数である。
hencky ひずみ $\epsilon $ は、変形勾配を用いて
\begin{align}
\epsilon = \frac{1}{2} \ln(FF^T)
\end{align}
と定義される。
具体的な応力を求めるため、
\begin{align}
B = FF^T
\end{align}
と置く。また、$B$ の転置を取ると
\begin{align}
B^T = (FF^T)^T=(F^T)^TF^T=FF^T=B
\end{align}
なので $B$ は対称テンソルである。 $B$ を変形勾配で微分すると
\begin{align}
\frac{\partial B_{ij}}{\partial F_{kl}} = \frac{\partial}{\partial F_{kl}} F_{im}F_{jm} = \delta_{ik} F_{jl} + \delta_{jk} F_{il}
\end{align}
となり、Piola-Kirchhoff 応力テンソルは、
\begin{align}
P_{kl} &= \frac{\partial \psi}{\partial F_{kl}}=\frac{\partial \psi}{\partial B_{ij}} \frac{\partial B_{ij}}{\partial F_{kl}} = 2 \frac{\partial \psi}{\partial B_{kj}}F_{jl} \\
P &= 2 \frac{\partial \psi}{\partial B}F
\end{align}
と書ける。
対数の行列の微分は
\begin{align}
\frac{\partial \ln B}{\partial B}B =B^{-1}B = I
\end{align}
と計算できるので、
\begin{align}
\frac{\partial \psi}{\partial B}B=\frac{1}{2}\frac{\partial \psi}{\partial \epsilon}\frac{\partial \ln B}{\partial B}B=\frac{1}{2}\frac{\partial \psi}{\partial \epsilon} = \frac{1}{2}D\colon\epsilon
\end{align}
である。Cauchy 応力は、
\begin{align}
\sigma = \frac{1}{J} P F^{T} \ \ , \ \ J = \mathrm{det}F
\end{align}
なので、 hencky ひずみにおけるCauchy 応力は、
\begin{align}
\sigma = \frac{1}{J} PF^T = \frac{2}{J} \frac{\partial \psi}{\partial B}FF^T= \frac{2}{J} \frac{\partial \psi}{\partial B}B = \frac{1}{J} D\colon\epsilon
\end{align}
と書くことができる。具体的には、
\begin{align}
J\sigma_{kl} = \lambda\delta_{kl}\epsilon_{mm} + 2\mu\epsilon_{kl}
\end{align}
と書けるので、これは、線形弾性体における応力である。
具体的な計算法 1 （Cauchy 応力を求める際に便利な方法）
変形勾配が求まり、hencky ひずみが計算できれば、線形弾性体のようにCauchy 応力を求めることができる。しかし、対数の中に行列 $B$ が入っており
\begin{align}
\epsilon = \frac{1}{2} \ln(FF^T) = \frac{1}{2} \ln(B)
\end{align}
の形なので単純には計算できない。しかし、行列 $B$ は対称行列なので、直交行列 $R$ を使い
\begin{align}
B = R^{-1}\Sigma R
\end{align}
と分解できる。$\Sigma$ は、$B$ の固有値が対角上に並んでいる対角行列である。よって、対数を $I$ の周りでテイラー展開を行い、べき展開したとき行列 $R^{-1},R$ が打ち消しあい
\begin{align}
\ln(B) &= \sum_{n=1} \frac{(-1)^{n+1}}{n} \left(B-I\right)^n \\
&= \sum_{n=1} \frac{(-1)^{n+1}}{n} \left(R^{-1}\Sigma R-R^{-1}R\right)^n \\
&= \sum_{n=1} \frac{(-1)^{n+1}}{n} R^{-1}\left(\Sigma-I\right)^n R \\
&= R^{-1}\ln(\Sigma) R
\end{align}
となる。$\Sigma$ は対角行列なので、$\ln(\Sigma) $ の計算は、対角成分に単純に対数を掛ければよい。
Cauchy 応力は、
\begin{align}
\sigma_{ij} &= \frac{1}{J} D_{ijkl}\epsilon_{kl} \\
&= \frac{1}{2J} D_{ijkl}\left(R^{-1}\ln(\Sigma) R\right)_{kl}
\end{align}
と計算できる。
具体的な計算法 2 （Piola-Kirchhoff 応力を求める際に便利な方法）
最後に、特異値分解によって計算する方法を紹介する。
Cauchy 応力は、行列 $B$ で記述できたので、行列 $B$ を対角化することで計算することができた。Piola-Kirchhoff 応力を求めたい場合、変形勾配を特異値分解した方が計算しやすい。
変形勾配を特異値分解すると、ユニタリー行列 $U,V$ を使い
\begin{align}
F = U\Sigma V^T
\end{align}
と分解できる。ユニタリー行列の特性を使えば、
\begin{align}
FF^T = U\Sigma V^T V \Sigma U^T = U\Sigma^2 U^T
\end{align}
と計算できるので、hencky ひずみが
\begin{align}
\epsilon &= \frac{1}{2}\ln(FF^T) = \frac{1}{2}\ln(U\Sigma^2 U^T) \\
&= \frac{1}{2}U\ln(\Sigma^2) U^T\\
&= U\ln(|\Sigma|) U^T
\end{align}
と求まる。また
\begin{align}
\frac{\partial \psi}{\partial \epsilon_{kl}} &= D_{klij}\epsilon_{ij}\\
& =\lambda\delta_{kl}\epsilon_{mm} + 2\mu\epsilon_{kl} \\
\frac{\partial \psi}{\partial \epsilon_{}} &=\lambda(\mathrm{Tr}\epsilon) I+ 2\mu\epsilon
\end{align}
であり、
\begin{align}
\frac{\partial \epsilon}{\partial F} &= \frac{1}{2} \frac{\partial \ln B}{\partial B } \frac{\partial B }{\partial F} = B^{-1}F \\
&=\left(U\Sigma^2 U^T\right)^{-1}U\Sigma V^T \\
&=U\Sigma^{-2} U^T U\Sigma V^T \\
&=U\Sigma^{-1} V^T
\end{align}
と計算できる。最終的に、Piola-Kirchhoff 応力は
\begin{align}
P &= \frac{\partial \psi}{\partial F} = \frac{\partial \psi}{\partial \epsilon}\frac{\partial \epsilon}{\partial F} \\
&= \left(\lambda(\mathrm{Tr}\epsilon) I+ 2\mu\epsilon \right)U\Sigma^{-1} V^T \\
&=\lambda\mathrm{Tr}\left(U\ln(|\Sigma|) U^T \right)U\Sigma^{-1} V^T +2\mu U\ln(|\Sigma|) U^T U\Sigma^{-1} V^T \\
&=U\left(\lambda\mathrm{Tr}\ln(|\Sigma|) \Sigma^{-1} + 2\mu \ln(|\Sigma|) \Sigma^{-1} \right)V^T
\end{align}
と計算できる。（途中でトレースの中の行列を巡回させて計算している。）
最後に
hencky モデルについて簡単にまとめたが、変形速度勾配 $L$ を求めて、
\begin{align}
L = \dot{F}F^{-1}
\end{align}
ひずみ速度
\begin{align}
\dot{\epsilon} = \frac{1}{2}\left(L+L^T\right)
\end{align}
を計算し、ひずみ（または応力）を更新していく方法が簡単かもしれない。（上手くいくかどうかは知らないが）