戻る
三次元移流拡散方程式の
有限要素解による各種変形要素の精度評価
（ 誤差解析手法によるアプローチ ）
( なお本文は下記からの一部抜粋および補遺論文である。）
（ 邵・松田・大河内, 三次元移流拡散方程式の有限要素解法における各種変形要素の精度評価
（誤差解析手法によるアプロ−チ), 日本機械学会論文集（B編）, Vol. 63-606 (1997) , pp. 444-451. )
１ 緒 言
非定常移流拡散方程式は熱流体解析の中で最も基本的な方程式の一つであり、
その効率的な数値解法の研究は重要である。著者らは二次元(1)(2)及び三次元(3)
移流拡散問題の有限要素解や差分解について、誤差解析手法により補正した
数値解法の提案を行ってきた。
有限要素法は本質的に汎用性をもつ数値解法であり、いろいろな要素が用いられている。
しかし各種の形状や変形要素の計算精度に関する理論的評価については、
二次元要素(5)に関しては著者らによる研究例あるものの三次元要素に関する研究例は
ほとんど無さそうである。この背景のもとに本論文は三次元移流拡散問題での
正四面体要素や、四面体及び六面体の不等長要素についての計算精度を理論的に
評価した。最後に数値実験により、誤差解析結果と数値計算結果との対応関係に
ついて検討した。 なお、本論文では次の三次元非定常移流拡散方程式を対象としている。
･･･(１)
ただし、 t : 時間，C : 濃度， u, v, w : x, y, z 方向の流速，
Kx, Ky, Kz : x, y, z 方向の拡散係数である。
２ 誤差解析手法による数値解法の補正方法
２．１ 一般解についての定式化
流速 u, v, w 及び拡散係数 Kx, Ky, Kz が定数である場合に式(１)は変数分離法によって
解くことができる。その一般解はフーリエ級数で表現され、その成分の一つは次式で表される。
･･･(２)
ただし、 A : 複素定数，βx,βy,βz : x, y, z 方向の波数である。
ここで、t=mDt，
, 
, 
，
, 
, 
，
, 
, 
とすれば、式(２)より次式を得る。
ただし m : 時間方向の添字，
, 
, 
: x, y, z 方向の添字，
, 
, 
: x, y, z 方向のメッシュ幅，Dt : タイムステップ である。
･･････････(３)
･････････････････････････････････････････････(４)
･････････････････････(５)
ここで 
は時刻 t=mDt での座標 , , における濃度であり、
は一般解における増幅因子である。
２．２ 数値解についての定式化と、その補正方法
一方、式(１)にある数値解法を適用し、変数 に 
を
代入して整理すると、数値解の増幅因子ξは一般に次式より表される。
･･･････････････････････････････････(６)
ここで、
：
の関数，
： 
と の関数，
： 
と の関数である。
ここで精度の向上のために、拡散項と移流項に補正係数ｆおよびｇを導入し、
次式のように補正したξを考える。
････････････････････････････････(７)
したがって、補正係数ｆとｇは式(４)と式(７)とを等値することにより求められるが、この際
に依存しない係数を求めるために、
の極限を考えると次式を得る。
･･････････････････････････････････(８)
ただし、
したがって、逆に式(８)で求めたｆとｇを式(７)に適用することによって、
数値解におけるξを補正することができる。
３ 有限要素法の陰解法による定式化とその補正
３．１ ガレルキン法による定式化
式(１)をガレルキン法により定式化すると、式(９)の関係式を考慮して
マトリックス形で式(10)のように表すことができる。
･･･････････････････････････････････(９)
･･････････････････････････････(10)
ただし、
，
θ：タイムスキームパラメータ，Ve：一要素領域，Ni：形状関数，M：総要素数である。
（１）正四面体一次要素の場合[ＧＲ(C)]
図１(a)に示した要素の場合、形状関数を四面体一次要素とすると
式(６)での 
は次式のように表される。ただし、
とする。
図１ 各種三次元要素
，
したがって、式(８)より補正係数として次式を得る。
･･････････････････････････････････････(11)
（２）不等長六面体一次要素の場合[ＧＨ(C)]
図１(b)に示す節点ＨＨの補正係数ｆ及びｇを示す。
･･････････････････････(12)
α＝β＝γ＝１のとき 
（３）不等長四面体一次要素の節点ＴＨ１の場合[ＧＴ１(C)]
図１(c)に示す節点ＴＨ１の補正係数ｆ及びｇを示す。
･･･････････(13)
α＝β＝γ＝１のとき、 
（４）不等長四面体一次要素の節点ＴＨ２の場合[ＧＴ２(C)]
図１(d)に示す節点ＴＨ２の補正係数ｆ及びｇを示す。
･･(14)
ただし、
α＝β＝γ＝１のとき 
････････････････････････････････(15)
，
である。
･････････････････････(16)
･････････････････････(17)
･･･････(18)
･･････(19)