工学系大学院生のブログ

2人の大学院生による雑記ブログ

第8-4回 一般座標系における熱伝導方程式の離散化[python, julia]



\begin{eqnarray}
\frac{\partial \hat{Q}}{\partial t} + \frac{\partial \hat{E}}{\partial \xi}+\frac{\partial \hat{F}}{\partial \eta}+\frac{\partial \hat{G}}{\partial \zeta}=0
\end{eqnarray}


第8-3回では,一般座標系における熱伝導方程式の離散化を行っていきます.

一般座標系の取り方は自由であるため,全ての格子において長さ1の直方体になる様に定義しています.(\(\Delta \xi =1, V=1\))





まず,デカルト座標系における熱伝導方程式は次の通りです.

$$ \frac{\partial }{\partial t}(\rho CT)-\frac{\partial}{\partial x} \left( \lambda \frac{\partial T}{\partial x} \right) – \frac{\partial}{\partial y} \left( \lambda \frac{\partial T}{\partial y} \right) – \frac{\partial}{\partial z} \left( \lambda \frac{\partial T}{\partial z} \right)=0 \ \ \ (1)$$


今回は一種類の材料のみを考えるので,熱拡散係数\(D\)を用いて次のようにまとめると,三次元における一般的な保存形の方程式は次のように書けます.

\begin{eqnarray}
\frac{\partial Q}{\partial t}+\frac{\partial E}{\partial x}+\frac{\partial F}{\partial y}+\frac{\partial G}{\partial z}=0 \ \ \ (2)
\end{eqnarray}


よって,保存形の方程式と比較すると次の結果が得られます.

\begin{eqnarray}
Q &=& T \ \ \ (3)\\
E &=& -\left( D \frac{\partial T}{\partial x} \right) \ \ \ (4)\\
F &=& -\left( D \frac{\partial T}{\partial y} \right) \ \ \ (5)\\
G &=& -\left( D \frac{\partial T}{\partial z} \right) \ \ \ (6)\\
D &=& \frac{\lambda }{\rho C}
\end{eqnarray}

第8-2回で求めた一般座標系におけるフラックスの評価方法について掘り下げていきます.


改めて,一般座標系におけるフラックスはデカルト座標系におけるフ
ラックスを用いて次のように書けます.

\begin{eqnarray}
\hat{E}_{i+1/2,j,k} = (A^\xi_x E+A^\xi_y F+ A^\xi_z G)_{i+1/2,j,k}\ \ \ (7)
\end{eqnarray}


境界における面積ベクトルは境界で定義されているため,計算することができます.


次に,デカルト座標系におけるフラックス(\(E_{i+1/2,j,k}\))についてチェインルールを用いて変形します.

\begin{eqnarray}
E = -D\frac{\partial T}{\partial x}\
= -D\left( \frac{\partial \xi}{\partial x}\frac{\partial T}{\partial \xi}+\frac{\partial \eta}{\partial x}\frac{\partial T}{\partial \eta}+\frac{\partial \zeta}{\partial x}\frac{\partial T}{\partial \zeta}\right) \\
= -D J \left( A^\xi_x\frac{\partial T}{\partial \xi}+A^\eta_x\frac{\partial T}{\partial \eta}+A^\zeta_x\frac{\partial T}{\partial \zeta}\right) \ \ \ (8)
\end{eqnarray}


まず,括弧の外側にある値は次のように評価します.ただし,今回の解析では,熱拡散係数は一定とするため\(D\)は定数です.また,ヤコビアンは体積比であるため,隣り合うセルの体積の平均を用いています.

\begin{eqnarray}
(D J)_{i+1/2,j,k} = D \left(\frac{2}{V_{i,j,k}+V_{i+1,j,k}} \right) \ \ \ (9)
\end{eqnarray}


括弧の中の第一項は次のようになります.

\begin{eqnarray}
A^\xi_x\frac{\partial T}{\partial \xi}_{i+1/2,j,k}=(A^\xi_x)_{i+1/2,j,k} (T_{i+1,j,k}-T_{i,j,k}) \ \ \ (10)
\end{eqnarray}


括弧の中の第二項は次のようになります.

\begin{eqnarray}
A^\eta_x\frac{\partial T}{\partial \eta}_{i+1/2,j,k}=\frac{(A^\eta_x)_{i+1/2,j,k}}{4} ((T_{i+1,j+1,k}-T_{i+1,j,k})+(T_{i+1,j,k}-T_{i+1,j-1,k})\\
+(T_{i,j+1,k}-T_{i,j,k})+(T_{i,j,k}-T_{i,j-1,k}))\\
= \frac{(A^\eta_x)_{i+1/2,j,k}}{4} ((T_{i+1,j+1,k}-T_{i+1,j-1,k})+(T_{i,j+1,k}-T_{i,j-1,k})) \ \ \ (11) \\
\end{eqnarray}


この式のイメージが少しでも分かりやすくなる様に,下の図を書きました.
現在は赤い部分のフラックスについて考えています.


この位置に置ける\(\eta\)方向の微分はつまるところ,青い矢印で出入りする部分の平均として処理しています.なので4で割っています.


同様に\((A^\eta_x)_{i+1/2,j,k}\)も青い矢印の部分の平均として計算します.

\begin{eqnarray}
(A^\eta_x)_{i+1/2,j,k} = \frac{1}{4}((A^\eta_x)_{i+1,j+1/2,k}+(A^\eta_x)_{i+1,j-1/2,k}+(A^\eta_x)_{i,j+1/2,k}+(A^\eta_x)_{i,j-1/2,k}) \ \ \ (12)
\end{eqnarray}


括弧内の第三項も同様です.

\begin{eqnarray}
A^\zeta_x\frac{\partial T}{\partial \zeta}_{i+1/2,j,k}=\frac{(A^\zeta_x)_{i+1/2,j,k}}{4} ((T_{i+1,j,k+1}-T_{i+1,j,k-1})+(T_{i,j,k+1}-T_{i,j-1,k-1})) \ \ \ (13) \\
(A^\zeta_x)_{i+1/2,j,k} = \frac{1}{4}((A^\zeta_x)_{i+1,j,k+1/2}+(A^\zeta_x)_{i+1,j,k-1/2}+(A^\zeta_x)_{i,j,k+1/2}+(A^\zeta_x)_{i,j,k-1/2}) \ \ \ (14)
\end{eqnarray}


式(9)~(14)の結果を式(8)に代入することで\(E_{i+1/2,j,k}\)を計算することが出来ます.


また,これらの考えを利用することで,\(F_{i+1/2,j,k}\),\(G_{i+1/2,j,k}\)も計算することができるため,式(7)から一般座標系におけるフラックスを求めることが出来ます.



離散化ができたので,次回は一般座標系における熱伝導方程式の結果についてお話します.

*コードは第8-5回にあります.

by hide

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です