長岡工業高等専門学校研究紀要

論第 52 巻（2016）

文

安定化気泡関数有限要素法による熱対流シミュレーション
今井 伸哉1・倉橋 貴彦2・片峯 英次3・河田 剛毅4
1 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
2 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
3 岐阜高専（Department of Mechanical Engineering, National Institute of Technology, Gifu College）
4 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Numerical Simulation for Thermal Convection
Based on Stabilized Bubble Function Finite Element Method
Shinya IMAI1, Takahiko KURAHASHI2, Eiji KATAMINE3 and Yoshitaka KAWADA4AbstractThe purpose of our study is to make a computer program to analyze thermal convection on the basis of
the stabilized finite element method using the bubble function element. As for the governing equation, theincompressible Navier-Stokes, the continuity, and the advection diffusion equations are employed, andthe stabilized finite element and the Euler methods are applied to discretize the governing equation in
space and time, respectively. The decoupled method is applied to divide the flow computation into thevelocity and the pressure fields. The bubble function element is applied to interpolate the velocity and thetemperature fields, and the linear triangular element is employed for the interpolation of the pressure field.
Key Words : Thermal convection, Incompressible Navier-Stokes equation, Boussinesq approximation,
Stabilized Finite Element Method, Bubble Function Element

１．序論
熱対流現象を解析することは，工学において重要である．例えば，工学において重要な対象の一つである熱機関は，熱エネルギーから作動流体を介して動力を得る．代表的な外燃機関の一つにスターリン
グエンジン（図-1）がある．この機関では，作動流体を，加熱・冷却により膨張・収縮させ駆動力を得る．そのため，このエンジンの作動流体を解析するには，熱と流れの連成を考慮する必要がある．

図-1 スターリングエンジンの実験用模型- 28 -今井 伸哉・倉橋 貴彦・片峯 英次・河田 剛毅

連成解析は求めるべき未知変数が多いため，計算コストを削減することは特に重要になる．有限要素法に基づいた解析において，計算コストおよび結果は，使用する要素タイプによる影響を受ける．本研
究では，気泡関数要素に着目する．二次元空間を三角形領域に分割するとき，気泡関数要素は，要素重心点に気泡関数点と呼ばれる節点を三角形一次要素に付加した 4 節点で各要素を補間する要素である．
一方，三角形二次要素は，要素内を 6 節点で補間する．ゆえに同じ要素数であれば，気泡関数要素は二次要素よりも計算コストが小さくなると考えられる．
気泡関数要素に関する研究について，安定化有限要ui x, t 0   uˆi ( x ) in  T x, t 0   Tˆ ( x ) in (4)on T  q
u x  u y  0on TT  Tˆ αT n  qˆ  0 on ,j jq(5)
ここで，t0は初期時刻， x = ( x,y ) は位置ベクトルを表す．なお，本解析において圧力は，領域全体に対する平均値が 0 になるように一意に定める．

素法との等価性については活発に議論されているが，他の要素タイプとの結果の比較はあまり例がなく，その例も定常の粘性流れ場に限られる

２.２ 支配方程式に対する気泡関数要素を用いた安．

定化有限要素法に基づく離散化1), 2), 3)

そこで本研究では，非定常熱対流問題に対して気気泡関数要素を用いた解析において，流れ場には泡関数要素を用いた安定化有限要素法に基づいて有流速場と圧力場を分離して解析する分離型解法を用

限要素方程式を導出し，それに基づく数値解析プロ

いる 4)．圧力場は，圧力ポアソン方程式を用いて解グラムによる解析結果を二次要素を用いた解析結果析を行う．圧力ポアソン方程式は，式(1)にEulerのと比較する．

前進差分法を適用した式(6)に対し発散を取り，式(2)に示す連続の式を考慮することにより得られる(式(7))．

２．支配方程式に対する離散化
uin+1 = uin  Δt (u nj ui,nj +２.１ 支配方程式
Γ = ΓT∪Γq で囲まれた二次元領域 Ω において，流 νui,nj + u nj,i , j + βg i T n  T0 )

れ場と温度場の連成場である自然対流場の解析を行う．本論文において用いる支配方程式を式(1)～(3) 1P,iin+1 =  ρ ui,ni  u njui,nj ,i  βgiT,in  Δtに示す．

u i  u j u i , j 1P,i  (u i , j  u j ,i ) , j  g i (T  T0 )  0 (1)ui ,i  0(2)T  u j T, j  T, jj  0(3)1 n+1P,iρ(6)(7)

なお，式(7) では粘性項に相当する項は，他の項より十分小さいとして無視する．
式 (7)，(6)，(3)に対し，支配方程式の空間方向の離散化を行う．各式に対し，重み関数 P*，ui*，T*を掛けて領域積分し，Greenの定理を用いることで，支配方程式の弱形式（式(8)～(10)）が得られる．

各式はそれぞれ，Boussinesq近似を適用した非圧縮
Navier-Stokes方程式，連続の式，エネルギー式である．ここで ui は流速 ( i = x,y )，P は圧力，T は温度，
T0 は基準温度，ρ は密度，ν は動粘性係数，β は熱膨張係数，gi は重力加速度，α は温度拡散係数を表す．なお，本論文の添え字表示では，偏微分表示(  ),i=∂(  )/∂xi および Einstein の総和規約を使用し，T などの変数上部の点はその変数の時間に対する偏導関数を表している．初期条件および境界条件は式(4)，(5)の様に定義する．
- 29 - P P dΩ =Ωe*,in+1,i

 1

ρ
P*ui,ni dΩ + P,i*u njui,nj dΩ  βg i P*T,in dΩ  (8)ΩeΩe Δt Ωe ρ  P u u n dΓ   P P n dΓ*Γenjni, j*iΓen+1,ii

安定化気泡関数有限要素法による熱対流シミュレーション



ここで，式(13)~(15)における係数行列は式(16)のよ u u dΩ =  u u dΩ  Δt  u u u dΩΩe*in+1iΩe

*ini



Ωe

*injni, j

うに書くことができる．

1ui,*i P n+1dΩ + ν ui,* j ui,nj dΩΩeρ Ωe(9)

[C ]  N N T d
,i,i
 ii

T
Si    N N B ,i  dDi u n = N ,i  N B T u j  N B, j T dΩ
Ωe

 nTn
B = ρu j N N B, j  ui ni dΓΓT+ N N , j  P n+1 ni dΓΓM   N B N B T dΩTTnn Au j  =  N B N B  u j  N B, j  dΩT H ji  N B, j N B,i  dΩt n = N t n dΓ
B i
 iΓ

qˆ= N B qˆdΓ
Γ





+ ν u u dΩ + βg i u T  T0 dΩ  u t dΓ Ωe
ΩeΓe




*i, j

nj,i

*i

n

*i i

 T TdΩ +  T u T dΩ + α T T dΩ +  T qˆdΩ = 0 (10)**ΩejΩe,jΩe*,j ,j*Γe

ここで，ti は表面力， q̂ は熱流束を表し，次式で与えられる．

 n  1 n+1nn ti =   P δij + νui, j + u j,i n j ρqˆ = αT n,j j    (11)



(16)




流速 ui，圧力 P，温度 T および，それぞれの変数に対する重み関数 ui*，P*，T* に対し，次に示す補間関数を適用する．

次に，時間微分項のある式(15)に対してθ法を適
ui = N1ui1 + N 2ui 2 + N 3ui 3 + N 4u~i 4 = N B  ui ~TT = N1T1 + N 2T2 + N 3T3 + N 4T4 = N B  T T P = N1 P1 + N 2 P2 + N 3 P3 = N  P *T****~*
ui = N1ui 1 + N 2ui 2 + N 3ui 3 + N 4ui 3 = N B  ui ~*
 *T*
**
*
T = N1T1 + N 2T2 + N 3T3 + N 4T4 = N B  T 
 P * = N P * + N P * + N P * = N T P * 
1 12 23 3

用し，時間方向に離散化をする．θ法では，状態変T数およびその時間微分を次のように近似する5)．
(12)

ここで，N1 ，N2，N3 ，N4 は形状関数であり，面積座標ηlを用いてηl=Nl (l=1,2,3)， N4=27η1η2η3と表せる．

  T n+1  T n T ΔtT   θ T n+1 + 1  θ T n 



(17)

式(17)を式(15)に代入することにより，時間方向にも離散化した有限要素方程式を得る（式(18)）．

この補間関数を支配方程式の弱形式に代入することM + θΔtAu + αH T (18)= M   1  θ Δt Au + αH T  qで，式(13)～(15)に示す有限要素方程式が得られる．

nn+1jjnnjj

C P =n+1ii

 1 ρ S i uin + Di u n uin  βg i S i T n  Δtn+ B 

次に数値的安定性のため，式(14)，(18)に，数値(13)粘性項を付加する．式(15)に対する数値粘性項は，次式で与えられる．

M u = M u n+1jnj 1T Δt Au u  S i  P n+1 + H jj uin ρnni

ν'(14) + H ji u nj + βg i M T n  T0   t in M T+ Au T + αH T +qˆ= 0jj(15)

N 4, j N 4, j dΩ =21  N 4 dΩ Ae τ SUPG  



(20)

ここで，Ae は要素面積である．式(20)の付加は，ν[Hjj]の4×4成分に式(19)の右辺項を付加することを意味する．ただし，τ SUPGはSUPG法に対する安定化パラメータであり，次式で与えられる4)．

- 30 -

今井 伸哉・倉橋 貴彦・片峯 英次・河田 剛毅

1  2  2  2 u  2  4ν  2  2 +  τ SUPG =    +   Δt   h   h 2   e   e  

(20)

ここで，Δt は時間増分量， u e は要素移流速度のノルム，he は要素長であり，本論文では次式を用いる．

u  + u e 2xe 2yhe  2 Ae0.6Point A

(21)

0.5

ue =(22)0.5

なお， uie は要素内の平均移流速度を表す．式(18)に対しては，式(19)，(20) の動粘性係数 ν を温度拡散0係数 α に書き換えればよい．なお，式(19)は定常移
流拡散場において，三角形気泡関数要素を用いた有図-2

解析モデル図限要素方程式を式変形によって気泡関数の自由度を削除した式と，三角形一次要素を用いたSUPG法による有限要素方程式が等価となるようにすることで導出される．

３．計算の流れ
本論文における，気泡関数要素を用いた安定化有限要素法による計算の流れを以下に示す．
図-3

Step1 メッシュ情報，計算条件，初期条件，境界条

有限要素メッシュ図

件の読み込む．
Step2 圧力ポアソン方程式(13)より，圧力 P(n+1)の解
表-1 物性値

析し，平均圧力を全節点から減じる（平均圧密度, ρ (kg/m3)
動粘度, ν (m2/s)
Prandtl 数, Pr (-)
基準温度, T0 (deg C)
熱拡散係数, α (m2/s)
体積膨張率, β (1/K)
重力加速度, gi (m/s2)

力を 0 とする）．
Step3 数値粘性項(式(19))を Navier-Stokes 方程式(14)に付加した式より，流速 ui(n+1)の解析する．
Step4 数値拡散項を移流拡散方程式(18)に付加した式より，温度 T(n+1)の解析する．
Step5 流速 ui(n+1)，温度 T(n+1)を更新して，Step2 に戻る．

表-2

節点数
要素数
時間増分量, Δt (s)時間ステップ数
θ 法の θ の値

４．数値解析例
４.１ 解析条件
本論文において対象とする数値解析のモデル図，

1000
1.0×10-6
7.01
0
ν/Pr0.000206
gx= 0, gy= - 9.8

計算条件

2601
5000
0.001
300,0000.5 (Crank Nicolson 法)

有限要素メッシュ図，各種物理量および計算条件を図-2，図-3，表-1，表-2に示す．

４.２ 解析結果
前節に示す条件下における，気泡関数要素を用い- 31 -

安定化気泡関数有限要素法による熱対流シミュレーション

た安定化有限要素法に基づく解析結果を，二次要素を用いた解析結果とともに以下に示す．

(b) 三角形二次要素
図-3 流速ベクトルと温度の分布図（時刻t = 300 s）(a)

三角形気泡関数要素(a)

三角形気泡関数要素(b) 三角形二次要素
図-2 流速ベクトルと温度の分布図（時刻t = 30 s）(b) 三角形二次要素
図-4 圧力の分布図（時刻t = 300 s）(a)

三角形気泡関数要素- 32 -

今井 伸哉・倉橋 貴彦・片峯 英次・河田 剛毅

数要素を用いた解析では，流れ場は分離型解法を用い，時間方向にはEulerの前進差分法，温度場は．
時間方向にはθ法を用い，空間方向にはいずれも安定化有限要素法に基づいて離散方程式を導出した．
二次要素（流れ場：直接法，Galerkin法, 完全陰解法，
温度場：Galerkin法，θ法）を用いた数値解析と解析結果を比較した．これらの比較より，定常状態において両者はおおよそ同様な分布を示すが，気泡関数要素を用いた解析の方が流速の大きさが小さく，過渡期において気泡関数要素を用いた解析の方が，温度の拡散が早いことが得られた．今後は，数学的に(a)どちらの傾向が正しいのか調査を行う予定である．

温度 Tの時間変化

謝辞：本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバ
SR16000を使用し計算を行ったものである．計算機の使用について，センター関係者の方へ謝意を表す．
参考文献
1)

松本 純一, 川原 睦人, 気泡関数要素を用いた GLS 型
安定化有限要素法による非圧縮性粘性流体解析, 第14 回数値流体力学シンポジウム, E06-4, 2000 (in CDROM).

(b) 流速の大きさ|u| = (ux2+uy2)1/2の時間変化
2)

図-5 座標 ( x, y ) = ( 0.6, 0.5 ) における状態量の時間変化

梅谷 信行, 非圧縮 Navier-Stokes 方程式に対する MINI
要素と安定化有限要素法の関係, 第 19 回数値流体力（P2: 二次要素，P1b: 気泡関数要素）

学シンポジウム, A4-1, 2005.

なお，二次要素を用いた解析において，流れ場と温3)

限要素法解析における要素の選択について, 土木学会

度場の連成は弱連成によって考慮し，流れ場は直接

年次学術講演会講演概要集, Vol.45, No.02, pp.568-569,

法および完全陰解法，温度場はθ 法を用い，空間方向の離散化には，Galerkin法をそれぞれ用いた．実

石毛 立也, 川原 睦人, 定常 Navier-Stokes 方程式の有1990.
4)

際の数値計算にはFreeFem++5)6)7)を用いた．

日本計算工学会流れの有限要素法研究委員会編，続・有限要素法による流れのシミュレーション, シュ

図-2，図-3 は各時刻における温度分布および流速分布を示し，図-4 は圧力分布を示す．図-3 ，プリンガー・ジャパン, 2008.
5)

図-4より，定常状態において両者は同様な分布を示
す こ と が わ か る ． ま た ，図 -5 は 図 -1 の Point A( x,y ) = ( 0.6, 0.5 ) における状態量の時間変化を示して

F. Hecht, New development in FreeFem++. Journal ofNumerical Mathematics, Vol.20, No.3-4, pp.251-266,
2012

6)

日本応用数理学会監修, 大塚厚二, 高石武史, 有限要素法で学ぶ現象と数理―FreeFem++数理思考プログラおり，気泡関数要素を用いた解析の方が，温度の拡ミング―, 共立出版, 2014

散が早く，流速が小さいことがわかる．
7)

片峯 英次, 今井 伸哉, 温度分布を規定する非定常自然対流場の形状同定問題の解法, 日本機械学会論文集 B
偏, Vol.82, No.833, 15-00578, 2016.

５．結論
本論文では，熱対流問題に対して気泡関数要素を用いた安定化有限要素法に基づいて離散化を行い，その有限要素方程式に基づく数値解析結果と二次要素を使用した数値解析結果の比較を行った．流れ場と温度場の連成は弱連成によって考慮した．気泡関- 33 -(2016. 8. 19 受付)