FEM for Heat Diffusion

熱拡散方程式の有限要素法による解法

Last Update:2009年05月17日

概要

熱が時間とともに，拡散していく方程式である熱拡散方程式の有限要素法解析の方法について述べる．まず熱拡散方程式の定式化を行い，弱形式に変換した後に，有限要素法の離散化を行う．軸対称問題における弱形式化と有限要素法離散化についても触れる．

有限要素法による定式化

支配方程式の導出

フーリエの法則より，熱流速ベクトル $\Phi$ は次のように温度場 T(x) の勾配によって表される．

$\Phi = -\lambda\nabla T$ 　　（１）

但し， $\lambda$ は熱伝導係数であり，温度勾配と熱流速を対応づける正定値対称テンソルである．

一方熱量保存の法則により，

$\rho c \frac{dT}{dt}= -\nabla\cdot\Phi + f$ 　　　　（２）

が成り立つ.但しは単位体積あたりの発熱量である．

(1)の両辺の発散をとり式(2)を代入する.

$\nabla\cdot\Phi = -\nabla\cdot(\lambda\nabla T)\\ \Leftright \rho c \frac{dT}{dt} = \nabla\cdot(\lambda\nabla T) + f$

となり，熱拡散方程式が得られる．

熱拡散方程式
$\rho c \frac{dT}{dt} = \nabla\cdot(\lambda\nabla T) + f$

$\lambda$ はテンソル量であるが，等方的な物性の場合は温度勾配と熱流速ベクトルの方向が一致するために，スカラー量として扱うことができる．また熱伝導率が一定であることを仮定すると次のようになる．

熱拡散方程式(熱伝導率が等方で一定)
$\rho c \frac{dT}{dt} = \lambda \nabla^2 T + f$ 　　　(3)

弱形式化

式(3)のように熱伝導係数が等方で一定であるとして，式展開を行う．

境界条件としては，次の２つを考える．

表面領域上でが与えられる
表面領域上で単位面積あたりの熱流入量が与えられる．

式(1)のフーリエの法則から

$Q=\lambda\frac{\partial T}{\partial n}$

となる．流出量ではなく流入量を境界条件として与えているのは，実際の問題では冷却よりも加熱する場合が多いという現実に即しているものである．

(3)の両辺に固定境界 S_1 上で０になる任意の仮想変分 $\delta T$ を掛けて解析領域全体で積分し，Gaussの発散定理を適応すると以下のようになる．

$\int_V \rho c \delta T\dot{T} dV= \int_V \delta T$\lambda \nabla^2 T + f $dV\;\;\;\;\forall\delta T \\ \Leftright \int_V \rho c \delta T \dot{T} dV= \int_V \nabla\cdot$ \lambda\delta T \nabla T$ - \lambda\nabla\delta T\nabla T+ \delta T f dV\;\;\;\;\forall\delta T\\ \Leftright \int_V \rho c \delta T\dot{T} + \lambda\nabla\delta T\nabla TdV= \int_V \delta T f dV + \int_S \b{n}\cdot$ \lambda\delta T \nabla T$ dS\;\;\;\;\forall\delta T\\ \Leftright \int_V \rho c \delta T\dot{T} + \lambda\nabla\delta T\nabla TdV= \int_V \delta T f dV + \int_S \delta T \lambda\frac{\partial T}{\partial n} dS\;\;\;\;\forall\delta T\\ \Leftright \int_V \rho c \delta T \dot{T} + \lambda\nabla\delta T\nabla TdV= \int_V \delta T f dV + \int_S \delta T Q dS\;\;\;\;\forall\delta T$

弱形式の熱伝導方程式
$\int_V\rho c \delta T\dot{T} + \lambda\nabla\delta T\nabla TdV= \int_V \delta T f dV + \int_S \delta T Q dS\;\;\;\;\forall\delta T$ 　　　(4)

有限要素法離散化

領域がメッシュ分割されているとしよう．式(4)の積分は各要素内での積分を全て足し合わせて得ることができる．つまり，

$\sum_e\int_{V_e} \rho c \delta T\dot{T} + \lambda\nabla\delta T\nabla TdV= \sum_e\int_{V_e} \delta T f dV + \sum_e\int_{S_e} \delta T Q dS\;\;\;\;\forall\delta T$ 　　　　(5)

ここで，ガラーキン法で変分を行うとして，要素内で次のように補間関数を用いてと $\delta T$ が離散化されているとする．

T = N^q T^q ， $\delta T = N^p \delta T^p$

また，領域の表面は次のように補間関数 $\bar{N}$ を用いて離散化されているとする.

$\delta T = \bar{N}^p \delta T^p$

これらを式(5)に代入すると以下のように，要素剛性行列，全体剛性行列を導出することができる．但し，式変形では要素内での番号付け $\{p,q\}\rightarrow \{a,b\}$ が成り立つとした．

$\sum_e\int_{V_e} \rho c N^p\delta T^pN^q\dot{T}^q + \lambda\nabla N^p\delta T^p\nabla N^qT^qdV= \sum_e\int_{V_e} N^p\delta T^p f dV + \sum_e\int_{S_e} \bar{N}^p\delta T^p Q dS\;\;\;\;\forall\delta T\\ \Leftright \sum_e$\delta T^p\int_{V_e} \rho c N^pN^q dV \dot{T}^q$ + \sum_e$\delta T^p\int_{V_e}\lambda\nabla N^p\nabla N^qdVT^q$= \sum_e$\delta T^p\int_{V_e} N^p f dV$ + \sum_e$\delta T^p \int_{S_e} \bar{N}^pQ dS$\;\;\;\;\forall\delta T\\ \Leftright\delta T^a$\sum_e M^{pq} $\dot{T}^b + \delta T^a$\sum_eA^{pq}$T^b= \delta T^a$\sum_e {F_{in}}^p_e $ + \delta T^a$\sum_e {F_{out}}^p_e $\;\;\;\;\forall\delta T\\ \Leftright\delta T^a M^{ab} \dot{T}^b + \delta T^a A^{ab}T^b= \delta T^a F_{in}^a + \delta T^a F_{out}^a \;\;\;\;\forall\delta T\\ \Leftright M^{ab} \dot{T}^b + A^{ab}T^b= F_{in}^a + F_{out}^a\\ \Leftright \[M\] \{\dot{T}\} + \[A\]\{T\}= \{F_{in}\} + \{F_{out}\}$

ここで， $M$ , $A$ , $\{F_{in}\}$ $\{F_{out}\}$ はそれぞれ全体質量行列，全体剛性行列,(体積力からくる)全体外力ベクトル,(表面力からくる)全体外力ベクトルである．

$\[M\]=M^{ab}=\sum_e M^{ab}_e$ ， $\[A\]=A^{ab}=\sum_e A^{ab}_e$ ， $\{F_{in}\}=F_{in}^a=\sum_e {F_{in}}^a_e$ ， $\{F_{out}\}=F_{out}^a=\sum_e {F_{ou}}^a_e$

これらは要素単位での次のような行列の足し合わせで得られることが分かる．

要素質量行列，剛性行列，外力ベクトル
$M_e^{pq}=\int_{V_e} \rho c N^pN^q dV$
$A_e^{pq}=\int_{V_e}\lambda\nabla N^p\nabla N^qdV$
${F_{in}}_e^p = \int_{V_e} N^p f dV$
${F_{out}}_e^p = \int_{S_e} \bar{N}^pQ dS$

得られた温度とその時間微分に関する連立一次方程式を陰的・陽的時間積分スキームを用いて解くことで，時間発展解が得られる．

軸対称問題

以下の図のように軸対称な３次元問題について考える，つまり，形状，初期条件，境界条件，物性値全てが軸対称な場合であり，解も軸対称となる．この場合は断面における２次元的な問題に帰着でき，３次元解析のコストを大幅に削減することができる．

弱形式化

軸対称な場合について，断面について２次元的に式(4)を離散化する．熱拡散方程式の弱形式は以下の通りであった．

$\int_V \rho c \delta T\dot{T} + \lambda\nabla\delta T\nabla TdV= \int_V \delta T f dV + \int_S \delta T Q dS\;\;\;\;\forall\delta T$

積分の変数を以下のように， $\{x,y,z\}\rightarrow \{r,\theta,z\}$ に変換しよう．

$\{\begin{array}{l}x\rightarrow rsin\theta\\ y\rightarrow rcos\theta\\ z\rightarrow z\end{array}$

この場合のヤコビアンの絶対値は

$\det J(r,\theta,z)=\det \[\begin{array} \sin\theta & r \cos\theta & 0 \\ \cos\theta & -r\sin\theta & 0 \\ 0 & 0 & 1 \end{array}\] = r$

となるので， $dxdydz = r drd\theta dz$ となる．また，この場合下の図のように，

$dS = rdsd \theta$

が成り立つので，これを上に代入すると，

$\int_\theta\int_r\int_z r\rho c \delta T\dot{T} + r\lambda\nabla\delta T\nabla T dzdrd\theta = \int_\theta\int_r\int_z \delta T f rdzdrd\theta + \int_\theta \int_s \delta T Q rdsd\theta\;\;\;\;\forall\delta T$

$\theta$ 方向に被積分関数は変化しないので， $\theta$ の積分を取り除くことができる．また，断面の２次元的な積分領域を,断面の境界をのように書くと，下のような軸対称問題についての弱形式化が得られる．

軸対称問題についての熱拡散方程式の弱形式
$\int_{v_e} r\rho c \delta T\dot{T} + r \lambda\nabla\delta T\nabla T dv = \int_{v_e} \delta T f rdv + \int_s \delta T Q rds\;\;\;\;\forall\delta T$ 　　　　(5)

有限要素法離散化

領域がメッシュ分割されているとしよう．式(4)の積分は各要素内での積分を全て足し合わせて得ることができる．つまり，

$\sum_e\int_{v_e} r\rho c \delta T\dot{T} + r\lambda\nabla\delta T\nabla T dv = \sum_e\int_{v_e} \delta T f rdv + \sum_e\int_{s_e} \delta T Q rds\;\;\;\;\forall\delta T$ 　　　　(5)

ここで，ガラーキン法で変分を行うとして，要素内で次のように補間関数を用いてと $\delta T$ が離散化されているとする．

T = N^q T^q ， $\delta T = N^p \delta T^p$

また，領域の表面は次のように補間関数 $\bar{N}$ を用いて離散化されているとする.

$\delta T = \bar{N}^p \delta T^p$

$\sum_e\int_{v_e} r\rho c N^p\delta T^pN^q\dot{T}^q + r\lambda\nabla N^p\delta T^p\nabla N^qT^qdv= \sum_e\int_{v_e} rN^p\delta T^p f dv + \sum_e\int_{s_e} r\bar{N}^p\delta T^p Q ds\;\;\;\;\forall\delta T\\ \Leftright \sum_e$\delta T^p\int_{v_e} r\rho c N^pN^q dV \dot{T}^q dv$ + \sum_e$\delta T^p\int_{v_e}r\lambda\nabla N^p\nabla N^qdVT^q dv$\\\;\;\;\;\;\;\;\;= \sum_e$\delta T^p\int_{v_e}r N^p f dv$ + \sum_e$\delta T^p \int_{s_e}r \bar{N}^pQ dS$\;\;\;\;\forall\delta T\\ \Leftright\delta T^a$\sum_e M^{pq} $\dot{T}^b + \delta T^a$\sum_eA^{pq}$T^b= \delta T^a$\sum_e {F_{in}}^p_e $ + \delta T^a$\sum_e {F_{out}}^p_e $\;\;\;\;\forall\delta T\\ \Leftright\delta T^a M^{ab} \dot{T}^b + \delta T^a A^{ab}T^b= \delta T^a F_{in}^a + \delta T^a F_{out}^a \;\;\;\;\forall\delta T\\ \Leftright M^{ab} \dot{T}^b + A^{ab}T^b= F_{in}^a + F_{out}^a\\ \Leftright \[M\] \{\dot{T}\} + \[A\]\{T\}= \{F_{in}\} + \{F_{out}\}$

$\[M\]=M^{ab}=\sum_e M^{ab}_e$ ， $\[A\]=A^{ab}=\sum_e A^{ab}_e$ ， $\{F_{in}\}=F_{in}^a=\sum_e {F_{in}}^a_e$ ， $\{F_{out}\}=F_{out}^a=\sum_e {F_{ou}}^a_e$

これらは要素単位での次のような行列の足し合わせで得られることが分かる．

要素質量行列，剛性行列，外力ベクトル
$M_e^{pq}=\int_{v_e} r\rho c N^pN^q dv$
$A_e^{pq}=\int_{v_e}r\lambda\nabla N^p\nabla N^qdv$
${F_{in}}_e^p = \int_{V_e} rN^p f dv$
${F_{out}}_e^p = \int_{S_e} r\bar{N}^pQ ds$

三角形一次要素を用いた場合の解析的な要素剛性行列の導出

三角形一次要素を用いて補間を行った場合，解析的に剛性行列を導出することができる．このような一次要素による線形な補間関数をの変わりにを使って表す．

軸までの距離も三角形内で一次的に変化するから

r = L^r r^r

というように，節点の値の補間によって得ることができる．これを用いて計算してみよう．

$\int_{v_e} (L^1)^a(L^2)^b(L^3)^c dv = \Delta\frac{2!a!b!c!}{(a+b+c+2)!}$ 　　　　(6);

が成り立つことを利用する．但し $\Delta$ は三角形の面積であるとする．

要素剛性行列

$A_e^{pq}=\int_\Omega L^rr^r\frac{\partial L^p}{\partial r}\frac{\partial L^q}{\partial r} + L^rr^r\frac{\partial L^p}{\partial z}\frac{\partial L^q}{\partial z} drdz\\\;\;\;\; = \frac{1}{3}(r^0+r^1+r^2)\{\frac{\partial L^p}{\partial r}\frac{\partial L^q}{\partial r} + \frac{\partial L^p}{\partial z}\frac{\partial L^q}{\partial z} \}\Delta$

結局この場合，重心位置での距離に２次元三角形のラプラス方程式の要素剛性を掛けたものになることが分かる．

要素質量行列

要素質量行列を計算する．式(6)から以下が成り立つ

$\int_{v_e} L^pL^qL^r dv \;\;\;= \{\begin{array}{ll}\Delta/10 & (p=q=r)\\ \Delta/60 & (p\ne q\ne r)\\ \Delta/30 & (p=q\ne r\;\;p\ne q=r,\;\;p=r\ne q) \end{array}$

これて次の要素質量行列は

$M_e^{pq}=\rho c $\int_{v_e} L^rL^pL^q dv$ r^r$

は次のように書ける．

$\{\begin{array}{l}M_e^{00} = \rho c \frac{\Delta}{60} (6r^0+2r^1+2r^2)\\ M_e^{11} = \rho c \frac{\Delta}{60} (2r^0+6r^1+2r^2)\\ M_e^{22} = \rho c \frac{\Delta}{60} (2r^0+2r^1+6r^2)\\ M_e^{01} = M_e^{10} = \rho c \frac{\Delta}{60} (2r^0+2r^1+1r^2)\\ M_e^{02} = M_e^{20} = \rho c \frac{\Delta}{60} (2r^0+1r^1+2r^2)\\ M_e^{12} = M_e^{21} = \rho c \frac{\Delta}{60} (1r^0+2r^1+2r^2)\end{array}$

円盤の冷却問題

「有限要素法による構造解析事例集」という本に載っている，円盤を周囲から冷やすした時の温度分布の時間履歴を調べるという例題を解いてみた．これを解くことでプログラムの精度検証を行った．問題の概要は，次のとおりである．

半径が1mで厚さが0.1mの円盤を初期状態t=0では全体が500℃として，t>0で円周を0℃として境界条件を設定する．円盤の物性値は鉄として，

熱伝導率(α)：48.0 W/(mK)
比熱(c) : 480 J/(kg℃)
密度(ρ) : 7.86*10^3  kg/m^3

のように設定する．このような熱拡散問題における，温度の時間履歴は厳密解が次のようになっている．

$T=T_0\sum_{n=1}^{\infty} A_n J_0$b_n\frac{r}{a}$e^{-p_n t}$

ただし，a:円盤の半径， T_0 :初期温度， J_0 :０次のベッセル関数，:０次のベッセル関数の根, J_1 :１次のベッセル関数， $A_n=\frac{2}{\beta_n J_1(\beta_n)}$ ， $p_n=\frac{\lambda}{c\rho}\frac{\beta_n^2}{a^2}$

さて，これを有限要素法で解析を行う．軸対称性を用いて下のように断面の２次元的な領域でのみ解析を行った．

補間関数は三角形１次要素として，時間積分は時間刻み dt=1 として，Crank-Nicolson法を用いた．ある時刻における断面の有限要素法解は次のとおり．

２時間後の温度分布の厳密解と有限要素法解の比較は次のとおり

有限要素法解が厳密解にぴったり合っていることが分かる

参考にしたもの

Books

有限要素法概説―理工学における基礎と応用 (FEM+BEM (3))	菊池文雄著
流れと熱伝導の有限要素法入門	矢川元基著
有限要素法による構造解析事例集　基礎問題から実務レベル問題まで	白鳥正樹，三好俊朗著

FEM for Heat Diffusion

目次

概要

有限要素法による定式化

支配方程式の導出

弱形式化

有限要素法離散化

軸対称問題

弱形式化

有限要素法離散化

三角形一次要素を用いた場合の解析的な要素剛性行列の導出

要素剛性行列

要素質量行列

円盤の冷却問題

参考にしたもの

Books

Links