定常伝熱解析定式化ガラーキン法による定式化 (軸対称問題)

解くべき方程式

ここでは，軸対称の定常の熱伝導方程式をがラーキン法を使い，有限要素法計算のために定式化します．解くべき偏微分方程式は，「定常伝熱解析:熱伝導方程式の説明」に示しています．定常状態の偏微分方程式は， \begin{align} \lambda\nabla^2 T+Q_V=0 \label{eq:SteadyThermal_Poisson} \end{align} です．ここで，\(T\) は温度 \(\unit{K}\)，\(Q_V\) は単位体積あたりの発熱量 \(\unit{J/(m^3\cdot sec)}\)，\(\lambda\)は熱伝導率 \(\unit{J/(m\cdot K\cdot sec)}\) です．

図1: 計算モデル

偏微分方程式からメッシュ分割

重み付き残差法

式(\ref{eq:SteadyThermal_Poisson})は，任意の関数 \(T^\ast\) で \begin{align} \int_V T^\ast\left(\lambda\nabla^2 T+Q_V\right)\diff V =0 \label{eq:Method_of_Weighted_Residuals} \end{align} が成り立つことと同じです．これは，重み付き残差法と呼ばれます．\(T^\ast\) が重み関数で，右辺が残差です．残差がゼロの時，元の微分方程式の解と一致します．この積分領域\(V\)は対象とする全ての領域です．後々便利なので，任意の関数\(T^\ast\)に， \begin{align} &\nabla T^\ast\cdot\vm{n} = 0& &\text{Vの境界} \label{eq:NeumannBoundaryCondition} \end{align} という制限を設けます．ここで，\(\vm{n}\)は境界領域で外側に向かう単位法線ベクトルです．このような制限を設けても，解が変わることはありません．この制限を設けグリーンの定理を用いると，式(\ref{eq:Method_of_Weighted_Residuals})は， \begin{align} \int_V \lambda\nabla T^\ast\cdot\nabla T\diff V + \int_S \lambda T^\ast\nabla T\cdot\vm{n}\diff S + \int_V T^\ast Q_V\diff V = 0 \end{align} と変形できます．

軸対称問題を解きますので，円柱座標系を使います．周方向 \(\theta\)は一定としますので，その微分はありません．すると， \begin{align} \int_V \lambda\left(\pdiffA{T^\ast}{r}\pdiffA{T}{r}+\pdiffA{T^\ast}{z}\pdiffA{T}{z}\right)\diff V + \int_S \lambda T^\ast\nabla T\cdot\vm{n}\diff S + \int_V T^\ast Q_V\diff V = 0 \label{eq:IntegralEquation_all} \end{align} が得られます．

領域のメッシュ分割

計算領域が複雑な場合，式(\ref{eq:IntegralEquation_all})の積分の計算は不可能でしょう．そこで，計算領域を図2に示すように単純な三角形に分割し，それぞれを積分することが考えられました．この操作は，しばしば メッシュ分割 と呼ばれます．そして，それぞれの三角形は エレメント と呼ばれます，領域$V$をエレメント\(v_i\)にメッシュ分割すると，式(\ref{eq:IntegralEquation_all})は， \begin{align} \sum_{i=1}^N\left[ \lambda\int_{v_i}\left(\pdiffA{T^\ast}{r}\pdiffA{T}{r}+\pdiffA{T^\ast}{z}\pdiffA{T}{z}\right)\diff V + \int_{s_i} T^\ast\lambda\nabla T\cdot\vm{n}\diff S + \int_{v_i} T^\ast Q_V\diff V\right] = 0\label{eq:IntegralEquation} \end{align} となります．エレメント内では熱伝導率 \(\lambda\) は一定なので，左辺第一項の積分の外側にすることができます．第二項も積分の外側にすることがでますが，後の計算のために積分の中に残し，勾配の前にします．エレメントは単純形状なので積分は容易になります．その一方で，積分の回数が増えます．単純計算(積分)を繰り返すことは，コンピューターが得意な作業です．

図2: 三角形エレメントを使ったメッシュ分割

エレメント内での温度の補間

ここでは，図2に示したある特定のエレメントを取り出し，そのエレメント内の温度を線形補間します．取り出したエレメントを図3に示します．エレメントが計算領域よりも十分小さい場合，もっと正確に言うと温度の変化が線形近似できるほどエレメントを小さくした場合，線形補間は良い近似になります．

図3: 三角形エレメント

この \(i\) 番目エレメント内の温度を線形近似関数， \begin{align} T^i(r,\,z)&=C_0^i + C_r^i r + C_z^i z\nonumber\\ &= \begin{bmatrix} 1 & r & z \end{bmatrix} \begin{bmatrix} C_0^i \\ C_r^i \\ C_z^i \end{bmatrix} \label{eq:linear_approximation} \end{align} で表します．これは，\((r,\,z,\ T)\)空間では，平面を表します．三角形エレメントの頂点では温度が決まっているので， \begin{align} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix}= \begin{bmatrix} 1 & r_1^i & z_1^i \\ 1 & r_2^i & z_2^i \\ 1 & r_3^i & z_3^i \end{bmatrix} \begin{bmatrix} C_0^i \\ C_r^i \\ C_z^i \end{bmatrix} \label{eq:LinearApproxmation_start} \end{align} を満たす必要があります．したがって，\((C_0^i,\,C_r^i,\,C_z^i)\)は， \begin{align} \begin{bmatrix} C_0^i \\ C_r^i \\ C_z^i \end{bmatrix}&= \begin{bmatrix} 1 & r_1^i & z_1^i \\ 1 & r_2^i & z_2^i \\ 1 & r_3^i & z_3^i \end{bmatrix}^{-1} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix}\nonumber\\ &=\cfrac{1}{\mathrm{det}(\vm{M})} \begin{bmatrix} r_2^i z_3^i-r_3^i z_2^i & r_3^i z_1^i-r_1^i z_3^i & r_1^i z_2^i-r_2^i z_1^i \\ z_2^i-z_3^i & z_3^i-z_1^i & z_1^i-z_2^i \\ r_3^i-r_2^i & r_1^i-r_3^i & r_2^i-r_1^i \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \end{align} となります．要するに，エレメント内の温度を線形近似することにしたので，三角形の頂点の温度から式(\ref{eq:linear_approximation})が決まります．ここで， \begin{align} \mathrm{det}(\vm{M})= \begin{vmatrix} 1 & r_1^i & z_1^i \\ 1 & r_2^i & z_2^i \\ 1 & r_3^i & z_3^i \end{vmatrix}=r_1^i z_2^i + r_2^i z_3^i + r_3^i z_1^i - r_3^i z_2^i -r_2^i z_1^i - r_1^i z_3^i \end{align} です．これは，行列式の性質から三角形エレメントの面積の二倍です．

これまでの結果をまとめると，\(i\)番目のエレメント内の線形近似した温度は， \begin{align} T^i(r,\,z) &=\cfrac{1}{\mathrm{det}(\vm{M})}\begin{bmatrix} 1 & r & z \end{bmatrix} \begin{bmatrix} r_2^i z_3^i -r_3^i z_2^i & r_3^i z_1^i-r_1^i z_3^i & r_1^i z_2^i -r_2^i z_1^i \\ z_2^i-z_3^i & z_3^i-z_1^i & z_1^i-z_2^i \\ r_3^i-r_2^i & r_1^i-r_3^i & r_2^i-r_1^i \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \end{align} となります．ちょっと式が長いので， \begin{align} T^i(r,\,z) &=\begin{bmatrix} N_1^i(r,\,z) & N_2^i(r,\,z) & N_3^i(r,\,z) \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \label{eq:LinearInterpolation_in_Element} \end{align} とします．\(N_m^i(r,\,z)\)は基底関数，あるいは形状関数と呼ばれ， \begin{align} N_1^i(r,\,z) & = \cfrac{1}{\mathrm{det}(\vm{M})}\left[(r_2^i z_3^i-r_3^i z_2^i)+(z_2^i-z_3^i)r+(r_3^i-r_2^i)z\right]\\ N_2^i(r,\,z) & = \cfrac{1}{\mathrm{det}(\vm{M})}\left[(r_3^i z_1^i-r_1^i z_3^i)+(z_3^i-z_1^i)r+(r_1^i-r_3^i)z\right]\\ N_2^i(r,\,z) & = \cfrac{1}{\mathrm{det}(\vm{M})}\left[(r_1^i z_2^i-r_2^i z_1^i)+(z_1^i-z_2^i)r+(r_2^i-r_1^i)z\right] \end{align} です．この形状関数には， \begin{align} &N_m^i(r_n^i,\,z_n^i)=\begin{cases} 0 & \text{\(m\neq n\)} \\ 1 & \text{\(m=n\)} \end{cases}& &\sum_{m=1}^3 N_m^i(r,\,z)=1& \end{align} という性質が有ります．

基底関数(形状関数)は，図2に示すメッシュ分割した時点で決定されます．そして，式(\ref{eq:IntegralEquation})を解くことで，各エレメントの頂点の温度がきまります．この温度が未知数です．元々，温度分布は連続的な量ですが，このように頂点の温度にすることで離散化されます．

エレメントの積分

重み関数

熱伝導の方程式を有限要素法で離散化するために計算を進めています．これまでに得られた結果は，次のとおりです．

1. 前々節では．熱伝導の偏微分方程式を積分方程式に変換し，エレメントに分割した．この節の結論は，式(\ref{eq:IntegralEquation})です．
2. 前節では，エレメント内の温度を線形補間しました．結論は，式(\ref{eq:LinearInterpolation_in_Element})です．

本節では，これらの結果を使いがラーキン法で離散方程式 (連立方程式) を求めます．

ガラーキン法では，重み関数\(T^\ast\)と偏微分方程式の解\(T\)は同じ基底関数を使います．すなわち \begin{align} T^i(r,\,z) &=\begin{bmatrix} N_1^i(r,\,z) & N_2^i(r,\,z) & N_3^i(r,\,z) \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \\ T^{\ast i}(r,\,z) &=\begin{bmatrix} N_1^i(r,\,z) & N_2^i(r,\,z) & N_3^i(r,\,z) \end{bmatrix} \begin{bmatrix} T_1^{\ast i} \\ T_2^{\ast i} \\ T_3^{\ast i} \end{bmatrix} \\ \end{align} です．この条件のもと，重み付き残差法の積分を行い，離散化を進めます．

積分(第一項:剛性マトリックス)

式(\ref{eq:IntegralEquation})の左辺の第一項の積分を計算します．そのために，\((\partial T^i/\partial r,\,\partial T^i/\partial z)\)を計算します．それは， \begin{align} \begin{bmatrix} \pdiffA{T^i}{r} \\ \pdiffA{T^i}{z} \end{bmatrix}&= \begin{bmatrix} \pdiffA{N_1^i}{r} & \pdiffA{N_2^i}{r} & \pdiffA{N_3^i}{r} \\ \pdiffA{N_1^i}{z} & \pdiffA{N_2^i}{z} & \pdiffA{N_3^i}{z} \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \\ &=\cfrac{1}{\mathrm{det}(\vm{M})} \begin{bmatrix} z_2^i-z_3^i & z_3^i-z_1^i & z_1^i-z_2^i \\ r_3^i-r_2^i & r_1^i-r_3^i & r_2^i-r_1^i \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} \end{align} となります．右辺は場所\(r,\,z\)の関数になっていないことに注意が必要です．エレメントの頂点の温度\((T_1^i,\,T_2^i,\,T_3^i)\)のみで決まり，一定の値です．このことは，線形補間の出発の式(\ref{eq:LinearApproxmation_start})から，容易に理解でいます．\((\partial T^{\ast i}/\partial r,\,\partial T^{\ast i}/\partial z)\)も同様に得られます．

この結果を利用すると，式(\ref{eq:IntegralEquation})の左辺の第一項の積分は， \begin{align} &\int_{v_i}\left(\pdiffA{T^\ast}{r}\pdiffA{T}{r}+\pdiffA{T^\ast}{z}\pdiffA{T}{z}\right)\diff V \\ &\qquad= \int_{v_i} \begin{bmatrix} \pdiffA{T^{\ast i}}{r} & \pdiffA{T^{\ast i}}{z} \end{bmatrix} \begin{bmatrix} \pdiffA{T^i}{r} \\ \pdiffA{T^i}{z} \end{bmatrix}\diff V \nonumber \\ &\qquad=\cfrac{1}{\left[\mathrm{det}(\vm{M})\right]^2}\nonumber\\ &\qquad\qquad\times\begin{bmatrix} T_1^{\ast i} & T_2^{\ast i} & T_3^{\ast i} \end{bmatrix} \begin{bmatrix} z_2^i-z_3^i & r_3^i-r_2^i \\ z_3^i-z_1^i & r_1^i-r_3^i \\ z_1^i-z_2^i & r_2^i-r_1^i \end{bmatrix} \begin{bmatrix} z_2^i-z_3^i & z_3^i-z_1^i & z_1^i-z_2^i \\ r_3^i-r_2^i & r_1^i-r_3^i & r_2^i-r_1^i \end{bmatrix} \begin{bmatrix} T_1^i \\ T_2^i \\ T_3^i \end{bmatrix}\nonumber\\ &\qquad\qquad\times\int_{v_i}\diff V\\ \end{align}

積分(第二項:熱の移動あるいは境界条件)

式(\ref{eq:IntegralEquation})の左辺の第二項の積分を計算します．この第二項は熱の流れを表します．温度の勾配に熱伝導係数を乗じた量は熱の移動(ベクトル量)を示します．これに境界面の法線方向の単位ベクトルとの内積をとり，面積分を行うと熱の移動の総和になります．すなわち \begin{align} \int_{s_i} \lambda\nabla T\cdot\vm{n}\diff S = \text{熱の流れ} \end{align} です．

図4に示すエレメントと境界条件を考えます．図中のエレメントの内，三角形の辺が隣り合う場合，熱の出入りは合計するとゼロになります．一方エレメントの辺から出て行く熱量は，その隣りの辺では同じだけの熱が入るからです．したがって， \begin{align} \sum_{i=1}^N \int_{s_i} T^\ast\lambda\nabla T\cdot\vm{n}\diff S&= \sum_{i=1} \int_{\Gamma_i} T^\ast\lambda\nabla T\cdot\vm{n}\diff S \label{eq:element_ThermalFlow} \end{align} となります．エレメントが隣接する部分の熱の流れは計算する必要がなく，問題で与えられている強化のみに着目して計算します．

伝熱問題の通常の境界条件は，図4に示すように(1)断熱あるいは(2)温度が与えられる — のふた通りが有ります．式で示すと， \begin{align} \nabla T\cdot\vm{n} &= 0& &\text{断熱の場合}& \nonumber\\ T&=T_0& &\text{温度が与えられる場合}& \end{align} です．前者(断熱)の場合，式(\ref{eq:NeumannBoundaryCondition})のために，自動的に成り立ちます (自然境界要件)．式(\ref{eq:element_ThermalFlow})の値をゼロにするだけなので，簡単です．

図4: 境界条件 (断熱 or 温度一定)

問題は，境界で温度が与えられる場合です．この状況で，式(\ref{eq:element_ThermalFlow})を計算します．もちろん，この面積分は境界の部分のみになります．図5に示す境界の $T=T_0$ を計算します． \begin{align} \int_{\Gamma_i} T^\ast\lambda\nabla T\cdot\vm{n}\diff S &= \lambda\begin{bmatrix} T_1^{\ast i} & T_2^{\ast i} & T_3^{\ast i} \end{bmatrix} \int_{\Gamma_i}\begin{bmatrix} N_1^i \\ N_2^i \\N_3^i \end{bmatrix} \nabla T\cdot\vm{n} \diff S\nonumber \\ \end{align}

図5: 温度一定の境界条件

積分(第三項:発熱)

ページ作成情報

参考資料

更新履歴

last update:2017/05/21 18:24:08