Posted at

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

\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$の座標をそれぞれ

\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)$は

\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$ の座標変化(変位)を

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

とする。

2点間の相対変位を

\delta \textbf{u} = \textbf{u}(r_2,t) - \textbf{u}(r_1,t) =

\begin{pmatrix}
\delta u \\
\delta v \\
\delta w
\end{pmatrix}

とすると、

\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}$は相対変位テンソルであり、

\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$で

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

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

\textbf{E} = \frac{1}{2} (\textbf{D} + \textbf{D}^T) \\

\Omega = \frac{1}{2} (\textbf{D} - \textbf{D}^T)

である。

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

参考文献:

1. 連続体の力学 (基礎物理学選書)

2. 弾性体の基礎


Juliaのコード

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

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

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はテンソル積を求める関数。

A =

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

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}

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点の座標を

\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点の座標を

\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}

とし、このときのひずみテンソルを求める。

グラフを描くと

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

青が時刻 t=0 での2点を、赤が時刻 t=1 での2点を表す。

青から赤への変化をあらわすひずみテンソルを求めていく。


計算結果

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}$の計算結果は

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$の計算結果は

Ω

3×3 Array{Float64,2}:

0.0 -0.0416667 0.125
0.0416667 0.0 0.166667
-0.125 -0.166667 0.0

となる。 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をインポートして使う。