Juliaで力学：ひずみテンソルを求める

```math
\textbf{D} = \textbf{E} + \Omega = 
\begin{pmatrix}
\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\
\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\
\frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z}
\end{pmatrix}
```

# 前置き

Juliaでひずみテンソルを計算してみました。
専門家ではないので、間違っているところがあるかもしれません。
何か気づかれるところがありましたら、コメントお願いします。

# 環境
MacOS Mojave
Julia v1.0.3

# 相対変位テンソルとひずみテンソル

3次元空間上に与えられた2点$\textbf{r}_1$と$\textbf{r}_2$の相対変位からひずみテンソルを求める問題。

時刻 t における$\textbf{r}_1$と$\textbf{r}_2$の座標をそれぞれ

```math
\textbf{r}_1(t)=
\begin{pmatrix}
x_1 \\
y_1 \\
z_1
\end{pmatrix}, \,
\textbf{r}_2(t)=
\begin{pmatrix}
x_2 \\
y_2 \\
z_2
\end{pmatrix}
```

とすると、2点間の距離$\textbf{r}(t)$は

```math
\delta \textbf{r} = \textbf{r}_2(t) - \textbf{r}_1(t) =
\begin{pmatrix}
\delta x \\
\delta y \\
\delta z
\end{pmatrix}
```

また、時刻 t と時刻 t+1 の間に起こった点 $\textbf{r}_1$ の座標変化（変位）を

```math
\textbf{u}(r_1,t) = \textbf{r}_1(t+1) - \textbf{r}_1(t)
```

とする。

2点間の相対変位を

```math
\delta \textbf{u} = \textbf{u}(r_2,t) - \textbf{u}(r_1,t) =
\begin{pmatrix}
\delta u \\
\delta v \\
\delta w
\end{pmatrix}
```
とすると、

```math
\delta \textbf{u} = 
\begin{pmatrix}
\delta u \\
\delta v \\
\delta w
\end{pmatrix}
= \textbf{D} \,
\begin{pmatrix}
\delta x \\
\delta y \\
\delta z
\end{pmatrix}
= \textbf{D} \, \delta \textbf{r}
```

と書ける。

ここで、$\textbf{D}$は相対変位テンソルであり、

```math
\textbf{D} =
\begin{pmatrix}
\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\
\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\
\frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z}
\end{pmatrix}
```

で求められる。

さらに、$\textbf{D}$は対称テンソル$\textbf{E}$と反対称テンソル$\Omega$で

```math
\textbf{D} = \textbf{E} + \Omega
```

と分解することができる。ただし、

```math
\textbf{E} = \frac{1}{2} (\textbf{D} + \textbf{D}^T) \\
\Omega = \frac{1}{2} (\textbf{D} - \textbf{D}^T)
```

である。

$\textbf{E}$がひずみテンソルであり、$\Omega$は回転を表すテンソルとなっている。

参考文献：
1. [連続体の力学 (基礎物理学選書) ](https://www.amazon.co.jp/%E9%80%A3%E7%B6%9A%E4%BD%93%E3%81%AE%E5%8A%9B%E5%AD%A6-%E5%9F%BA%E7%A4%8E%E7%89%A9%E7%90%86%E5%AD%A6%E9%81%B8%E6%9B%B8-%E4%BD%90%E9%87%8E-%E7%90%86/dp/4785321377)
2. [弾性体の基礎](http://www.hapis.k.u-tokyo.ac.jp/public/makino/materials/20081021_elastic.pdf)

# Juliaのコード

与えられた座標データからひずみテンソルを求める関数を実装。

$\partial$のところを$\delta$として考えて個々の要素を計算した。

```julia
function strain_tensor(r1_t1,r1_t2,r2_t1,r2_t2)
    δr = r2_t1 - r1_t1
    r1_δt = r1_t2 - r1_t1
    r2_δt = r2_t2 - r2_t1
    δu = r2_δt - r1_δt
    D = kron(δu,1/(δr)) #相対変位テンソル
    E = 1/2*(D+transpose(D)) #ひずみテンソル
    Ω = 1/2*(D-transpose(D)) #回転のテンソル
    return D,E,Ω
end
```

点$\textbf{r}_1$と$\textbf{r}_2$の時刻 t1 と時刻 t2 のそれぞれの座標を引数として、相対変位テンソル、ひずみテンソル、回転のテンソルを返す。

途中で使っている`kron`はテンソル積を求める関数。


```math
A =
\begin{pmatrix}
x_a \\
y_a \\
z_a
\end{pmatrix}, \,
B =
\begin{pmatrix}
x_b & y_b & z_b
\end{pmatrix}
```

```math
A \otimes B =
\begin{pmatrix}
x_ax_b & x_ay_b & x_az_b \\
y_ax_b & y_ay_b & y_az_b\\
z_ax_b & z_ay_b & z_az_b
\end{pmatrix}
```

```julia
A = [1.,2.,3.]
B = transpose([1.,4.,9.])
kron(A,B)
```
```
3×3 Array{Float64,2}:
 1.0   4.0   9.0
 2.0   8.0  18.0
 3.0  12.0  27.0
```

# 例題

時刻 t=0 での2点の座標を

```math
\textbf{r}_1(0) = 
\begin{pmatrix}
0.0 \\
0.0 \\
0.0
\end{pmatrix}, \,
\textbf{r}_2(0) = 
\begin{pmatrix}
0.2 \\
0.2 \\
0.2
\end{pmatrix}
```
時刻 t=1 での2点の座標を

```math
\textbf{r}_1(1) = 
\begin{pmatrix}
0.0 \\
0.0 \\
0.0
\end{pmatrix}, \,
\textbf{r}_2(1) = 
\begin{pmatrix}
0.3 \\
0.35 \\
0.15
\end{pmatrix}
```
とし、このときのひずみテンソルを求める。

グラフを描くと

```julia
using Plots
gr()

r1_t1 = [0,0,0]
r1_t2 = [0,0,0]
r2_t1 = [0.2,0.2,0.2]
r2_t2 = [0.3,0.35,0.15]

xyz_t1 = [r1_t1 r2_t1]
xyz_t2 = [r1_t2 r2_t2]
plot(xyz_t1[1,:],xyz_t1[2,:],xyz_t1[3,:],lw=3)
plot!(xyz_t2[1,:],xyz_t2[2,:],xyz_t2[3,:],legend=:none,
    xlims=(0.0,0.35),ylims=(0.0,0.35),zlims=(0.0,0.35),lw=3)
scatter!(xyz_t1[1,:],xyz_t1[2,:],xyz_t1[3,:],color="blue")
scatter!(xyz_t2[1,:],xyz_t2[2,:],xyz_t2[3,:],color="red")
```

![20190214_graph.png](https://qiita-image-store.s3.amazonaws.com/0/204284/41b090c5-825e-81fc-5acf-efed1b038496.png)

青が時刻 t=0 での2点を、赤が時刻 t=1 での2点を表す。
青から赤への変化をあらわすひずみテンソルを求めていく。

# 計算結果

```julia
function strain_tensor(r1_t1,r1_t2,r2_t1,r2_t2)
    δr = r2_t1 - r1_t1
    r1_δt = r1_t2 - r1_t1
    r2_δt = r2_t2 - r2_t1
    δu = r2_δt - r1_δt
    D = kron(δu,1/(δr)) #相対変位テンソル
    E = 1/2*(D+transpose(D)) #ひずみテンソル
    Ω = 1/2*(D-transpose(D)) #回転のテンソル
    return D,E,Ω
end

D,E,Ω = strain_tensor(r1_t1,r1_t2,r2_t1,r2_t2);
```

ひずみテンソル$\textbf{E}$の計算結果は

```julia
E
```
```
3×3 Array{Float64,2}:
 0.166667   0.208333    0.0416667
 0.208333   0.25        0.0833333
 0.0416667  0.0833333  -0.0833333
```

回転のテンソル$\Omega$の計算結果は

```julia
Ω
```
```
3×3 Array{Float64,2}:
  0.0        -0.0416667  0.125   
  0.0416667   0.0        0.166667
 -0.125      -0.166667   0.0 
```

となる。 Juliaは行列の出力が見やすくて良い感じ。


ひずみテンソル等から求められる各種パラメータも計算してみる。

```julia
using LinearAlgebra

print("体積膨張率は",round(tr(E),sigdigits=3),"\n\n")
print("x-y平面での回転は",round(Ω[2,1],sigdigits=3),"\n")
print("z-x平面での回転は",round(Ω[1,3],sigdigits=3),"\n")
print("y-z平面での回転は",round(Ω[3,2],sigdigits=3),"\n\n")
print("x-y平面でのひずみは",round(E[1,2],sigdigits=3),"\n")
print("z-x平面でのひずみは",round(E[1,3],sigdigits=3),"\n")
print("y-z平面でのひずみは",round(E[2,3],sigdigits=3),"\n")
```
```
体積膨張率は0.333

x-y平面での回転は0.0417
z-x平面での回転は0.125
y-z平面での回転は-0.167

x-y平面でのひずみは0.208
z-x平面でのひずみは0.0417
y-z平面でのひずみは0.0833
```

ちなみに、`tr`は$Trace$（トレース）。Juliaの標準の関数ではないので注意。
`LinearAlgebra`をインポートして使う。