論 文

三次元仮想領域FEMによる埋め込み物体を考慮した温度場の解析
倉橋 貴彦1・Askhar AIBOLAT1・河田 剛毅2
1 長岡技術科学大学大学院

機械創造工学専攻（Department of Mechanical Engineering, Graduate School of Nagaoka University of Technology）
2 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Numerical Analysis of Temperature Field Considering Embedded Body
Based on the Fictitious Domain FEM in Three Dimensions
Takahiko KURAHASHI1, Askhar AIBOLAT1 and Yoshitaka KAWADA2
Abstract
In this study, we present three dimensional temperature field analysis considering an embedded body
based on the fictitious domain FEM. As the governing equation, the steady and the unsteady heat transfer
equations are introduced, and the finite element method is applied to discretize the governing equations inspace. The Lagrange multiplier method is applied to consider boundary condition on the boundary of theembedded domain. In addition, the search of nodal points of embedded domain in the background meshes
is carried out in the finite element analysis, and this is one of the feature in this study.
Key Words : fictitious domain finite element method, temperature field, three-dimensional analysis

めの方法として，重合格子法1), 2)や埋め込み境界法3)，１．はじめにまた仮想領域法 4), 5)という方法がある．本論文では，を考慮した解析を実施することは多々あり，一般的仮想領域法を有限要素法(FEM)に適用した仮想領域FEM6)を使用し，重合メッシュにおける節点探査をな数値計算法を用いる場合，解析領域内における物三次元モデルに拡張することにより，埋め込み物体体をくり抜いた領域に対してメッシュ分割を行い，を考慮した三次元温度場の解析について検証を行う．

数値解析において，解析対象となる領域内に物体流れ場や温度場，応力場等に関する数値解析を行う．
しかしながら，物体の位置を動かして再度解析を実２．三次元仮想領域FEMによる定式化施する場合等，解析領域内の物体の位置を動かしながらその都度メッシュやグリッドを生成し，解析を実施する必要がある．このような問題を解決するた定式化の説明に際し，三次元モデルにおける定常1温度場の支配方程式（式(1)）を導入する．

𝑘𝑘 �

𝜕𝜕2 𝜙𝜙𝜕𝜕𝑥𝑥 2+

𝜕𝜕2 𝜙𝜙𝜕𝜕𝑦𝑦 2+

𝜕𝜕2 𝜙𝜙𝜕𝜕𝑧𝑧 2�=0𝜕𝜕𝜕𝜕𝜕𝜕𝜙𝜙𝑘𝑘(1)

𝜙𝜙4 は「式(3)の解（𝜙𝜙1 ～𝜙𝜙4 ）」になることがわかる．

辺に重み関数𝑤𝑤をかけ，要素領域Ω𝑒𝑒 において積分す

ここで，四面体一次要素内のある点の位置𝒙𝒙において温度𝜙𝜙が規定されており𝜙𝜙(𝒙𝒙) = 𝑓𝑓とする．形状関ることにより，式(2)に示す重み関数が得られる．
𝜕𝜕2 𝜙𝜙𝑒𝑒

𝜕𝜕2 𝜙𝜙

𝜕𝜕2 𝜙𝜙(7)式(7)は式(3)と等価であり，式(6)を最小とする𝜙𝜙1 ～ここに，kは熱拡散率，𝜙𝜙は温度を示す．式(1)の両∫Ω 𝑤𝑤𝑤𝑤 � 𝜕𝜕𝑥𝑥 2 + 𝜕𝜕𝑦𝑦2 + 𝜕𝜕𝑧𝑧 2 � 𝑑𝑑Ω = 0= ∑4𝑗𝑗=1�𝐻𝐻𝑘𝑘𝑘𝑘 𝜙𝜙𝑗𝑗 − 𝑞𝑞𝑘𝑘 � = 0数𝑁𝑁𝑖𝑖 を用いることで，式(8)のように記載することができる．

(2)

𝑁𝑁1 (𝒙𝒙)𝜙𝜙1 + 𝑁𝑁2 (𝒙𝒙)𝜙𝜙2 + 𝑁𝑁3 (𝒙𝒙)𝜙𝜙3 + 𝑁𝑁4 (𝒙𝒙)𝜙𝜙4 = 𝑓𝑓

式(2)に対してグリーンの定理を適用し，四面体要素による一次の補間関数を重み関数𝑤𝑤および温度𝜙𝜙に適用すると，式(3)の形の代数方程式が得られる．

(8)

今，式(8)を左辺に移項した{𝑵𝑵}𝑇𝑇 {𝝓𝝓} − 𝑓𝑓 = 0を，式(6)の最小化問題に対する制約条件として考える．
ラグランジュの未定乗数法を用いることで，式(9)

𝐻𝐻11
𝐻𝐻21
�
𝐻𝐻31
𝐻𝐻41

𝐻𝐻12
𝐻𝐻22
𝐻𝐻32
𝐻𝐻42

𝐻𝐻14 𝜙𝜙1𝑞𝑞1
𝐻𝐻24 𝜙𝜙2𝑞𝑞2� � � = �𝑞𝑞 �
𝐻𝐻34 𝜙𝜙33𝑞𝑞4
𝐻𝐻44 𝜙𝜙4

𝐻𝐻13
𝐻𝐻23
𝐻𝐻33
𝐻𝐻43のように書くことができる．式(9)において𝜆𝜆はラグランジュ未定乗数であり，式(9)はラグランジュ関(3)数と呼ばれる．
1𝐽𝐽∗ = {𝝓𝝓}𝑇𝑇 [𝑯𝑯]{𝝓𝝓} − {𝝓𝝓}𝑇𝑇 {𝒒𝒒} + 𝜆𝜆({𝑵𝑵}𝑇𝑇 {𝝓𝝓} − 𝑓𝑓)

ここに𝐻𝐻𝑖𝑖𝑖𝑖 (i=1～4，j=1～4)は有限要素法による定式2化により得られるものであり，式(4)により計算す(9)

ラグランジュ関数J*を最小とする𝜙𝜙1 ～𝜙𝜙4 およびラグることができる．

ランジュ未定乗数𝜆𝜆を求めるために各変数により微𝐻𝐻𝑖𝑖𝑖𝑖 = ∫Ω𝜕𝜕𝑁𝑁𝑖𝑖 𝜕𝜕𝑁𝑁𝑗𝑗

𝑒𝑒 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕+𝜕𝜕𝑁𝑁𝑖𝑖 𝜕𝜕𝑁𝑁𝑗𝑗
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕+𝜕𝜕𝑁𝑁𝑖𝑖 𝜕𝜕𝑁𝑁𝑗𝑗
𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕𝑑𝑑Ω

分すると式(10)が得られる．

(4)

𝜕𝜕𝐽𝐽 ∗�[𝑯𝑯]
𝜕𝜕𝝓𝝓� 𝜕𝜕𝐽𝐽∗ � = �{𝑵𝑵}𝑇𝑇�

式(4)において，𝑁𝑁𝑖𝑖 は四面体一次要素における形状関数を示しており，式(5)により計算することができる．

𝜕𝜕𝜕𝜕{𝑵𝑵} {𝝓𝝓}{𝒒𝒒}{𝟎𝟎}���−� �=� �
𝑓𝑓00
𝜆𝜆(10)

式(10)より結果として式(11)が得られ，式(11)を解く(5)

𝑁𝑁𝑖𝑖 (𝑥𝑥, 𝑦𝑦, 𝑧𝑧) = 𝑎𝑎𝑖𝑖 + 𝑏𝑏𝑖𝑖 𝑥𝑥 + 𝑐𝑐𝑖𝑖 𝑦𝑦 + 𝑑𝑑𝑖𝑖 𝑧𝑧

ことにより，四面体要素内のある点の位置𝒙𝒙における温度𝜙𝜙(𝒙𝒙) = 𝑓𝑓を満たすように温度𝜙𝜙1 ～𝜙𝜙4 を求め𝑎𝑎𝑖𝑖 ，𝑏𝑏𝑖𝑖 ，𝑐𝑐𝑖𝑖 ，𝑑𝑑𝑖𝑖 は四面体要素の節点の各座標値によることができる．式(11)が三次元仮想領域FEMにお7)ける有限要素方程式である．

り計算されるものであり，式の詳細は文献 を参照して頂くことにする．今，式(3)を[𝑯𝑯]{𝝓𝝓} = {𝒒𝒒}と表示することにする．大括弧は行列，中括弧はベクト�ルを示している．この式に対して，式(6)のように関数Jを定義する．

[𝑯𝑯]{𝑵𝑵}𝑇𝑇{𝑵𝑵} {𝝓𝝓}{𝒒𝒒}���=� �𝑓𝑓0𝜆𝜆(11)非定常モデルの場合は，式(12)のように書くことが1

𝐽𝐽 = {𝝓𝝓}𝑇𝑇 [𝑯𝑯]{𝝓𝝓} − {𝝓𝝓}𝑇𝑇 {𝒒𝒒}2できる．時間微分については，差分法を適用するこ(6)とにより，離散化を行う．本論文の数値解析例では，クランク・ニコルソン法を適用し，数値実験を行う．

ここで，関数Jを𝜙𝜙1 ～𝜙𝜙4 により微分し，行列[𝑯𝑯]の対称性も考慮すると式(7)が得られる．式(7)において𝑘𝑘は1～4のそれぞれの値を示す．
2

𝜕𝜕{𝝓𝝓}[𝑯𝑯]� 𝜕𝜕𝜕𝜕 � + �
𝑇𝑇{𝑵𝑵}
0{𝑵𝑵} {𝝓𝝓}{𝒒𝒒}���=� �𝑓𝑓0𝜆𝜆(12)式(11)および式(12)の計算は通常の有限要素解析のプロセスと基本的には同様であるが，どこの外部領域の要素内に内部領域の節点が存在するか判定を行う必要がある．判定方法の詳細は次章にて述べる．

３．内部領域の節点の内外判定について
仮想領域FEMでは，式(8)に示す補間関数を用いて内部領域の節点の状態量を外部領域の要素の節点

外部領域の要素 ijkl と内部要素の P 点の一例

図-1に与える処理を行う．そのため，内部領域のある節点が，外部領域において対象とした要素の内側にあ４．数値解析例

るか，外側にあるかについての判定を行う必要がある．ここで，図-1に示すように外部領域の要素ijkl内に，内部領域の節点Pがある場合を考える．この本検討では，図-2に示す三次元領域に対して，ま場合，外部領域の四面体の平面全てにおいて式(13)ず定常温度場の問題に対してプログラムの妥当性にに示す不等式が全て正になる必要がある8)．

ついて検証を行う．外部領域の正六面体の構造の1辺のサイズは，4m，内部領域の正六面体の構造の1𝐴𝐴𝑥𝑥𝑝𝑝 + 𝐵𝐵𝑦𝑦𝑝𝑝 + 𝐶𝐶𝑧𝑧𝑝𝑝 + 𝐷𝐷 ≥ 0(13)

辺のサイズは，1mとし，各正六面体の中心の位置を原点とする．境界条件としては，外部領域の外側境界に0℃，内部領域の各節点に100℃を与える．そここに，式(13)の係数A～Dは式(14)～式(17)により計算するこ とができる．式 (14) ～式(17)におけ るの他の計算条件を表-1に示す．時間増分量，時間ス(i,j,k)は，(k,j,l)，(i,l.j)，(k,l,i)と変え，四面体の境界テップの条件は，非定常モデルに対する検討の際の面の四面全てにおいて式(13)の不等式を満たす場合条件である．非定常モデルにおける解析においては，に，対象とした要素ijklの内部に内部領域の節点Pが外部領域の全節点に温度零℃の初期条件を与える．

存在するということになる．

z𝑦𝑦𝑖𝑖
𝐴𝐴 = � 𝑦𝑦𝑗𝑗
𝑦𝑦𝑘𝑘𝑧𝑧𝑖𝑖
𝑧𝑧𝑗𝑗
𝑧𝑧𝑘𝑘

11�
1(14)𝑧𝑧𝑖𝑖
𝐵𝐵 = � 𝑧𝑧𝑗𝑗
𝑧𝑧𝑘𝑘𝑥𝑥𝑖𝑖𝑥𝑥𝑗𝑗
𝑥𝑥𝑘𝑘

11�1(15)𝑥𝑥𝑖𝑖
𝐶𝐶 = � 𝑥𝑥𝑗𝑗
𝑥𝑥𝑘𝑘

𝑦𝑦𝑖𝑖
𝑦𝑦𝑗𝑗
𝑦𝑦𝑘𝑘

11�1(16)

𝐷𝐷 = −𝐴𝐴𝑥𝑥𝑖𝑖 − 𝐵𝐵𝑦𝑦𝑖𝑖 − 𝐶𝐶𝑧𝑧𝑖𝑖

yx(17)

図-2

計算モデル図（外側：外部領域，内側：内部領域）3

表-1 計算条件

100

98,30422,065
1536
369
0.1
1.0
40

90
80
7060
Temperature[℃]

外部領域の総節点数
外部領域の総要素数
内部領域の総節点数
内部領域の総要素数
時間増分量Δt，sec.
熱拡散率k，m2/ sec.
時間ステップ50z=0のx-y断面およびx=0のy-z断面における温度の40分布図を図-3に示す．結果より，温度分布が適切に30求まっていることを確認できる．図-4はy=0，z=0におけるx軸上の温度分布を示しており，結果より，20左端において完全には100℃になっていないことが10わかる．これは，補間関数（式(8)）により境界条件を規定していることに起因しており，直接的に外0部領域の節点に境界条件の値を与えていないことを

11.25
1.5
1.75
Distance from origin[m]（ｙ＝０、ｚ＝０）意味している．
図-4

また，非定常モデルにおいてx=0のy-z断面におけ2y=0，z=0 における x 軸上の温度分布（定常モデルの結果）る温度の分布図を図-5に示す．10ステップ（t=1sec.）
おきに表示したものである．図-4に示す定常モデル
100の結果に時間進展とともに近づいている様子を確認
で き る ． ま た ， 図 -6 は 各 熱 拡 散 率 (k=1.0 ， 2.0 ，903.0m2/sec.)の条件において検討をしたものであり，Temperature[℃]設定した場合，他の熱拡散率の条件の結果に比べて温度が早めに拡散している様子を確認できる．以上の結果より，定常・非定常モデルの両モデルにおいて三次元の温度場解析に仮想領域FEMを適用し，適切に解を得られていることを確認できる．

60

Steady-state Heat transferequation

50
402010[℃]

図-34[s]30z0x3[s]70

分布を示している．結果として，熱拡散率を大きく

x2[s]

80

時刻t=1.0sec.においてy=0，z=0におけるx軸上の温度

y1[s]

y1

1.25
1.5
1.75
Distance from origin[m]（ｙ＝０、ｚ＝０）

図-5

y=0，z=0 における x 軸上の温度分布（非定常モデルの結果）

z=0 の x-y 断面および x=0 の y-z 断面における温度の分布図

423) Kino,C., Watanabe,T., Nishida,A. and Takemiya,H., Flow

100analysis around an in-line forced oscillating circular cylinder

k = 1[m²/s]

90using IB-method, Transactions of the Japan Society ofk = 2[m²/s]

80

Mechanical Engineers Series B, Vol.78, No.796, pp.60-73,
2012.

k = 3[m²/s]704) Glowinski,R., Pan,T.W., Hesla,T.I., Joseph,D.D., andPeriaux,J., A Fictitious Domain Method with Distributed

Temperature[℃]60

Lagrange Multipliers for the Numerical Simulation of50

Particulate Flow, Contemp. Math., Vol.218, pp.121-137,

401998.
5) Glowinski,R., Pan,T.W., Hesla,T.I., Joseph,D.D. and30

Periaux,J., A Fictitious Domain Approach to the Direct20

Numerical Simulation of Incompressible Viscous Flow pastMoving Rigid Bodies: Application to Particulate Flow,

100

図-6

J. comput. phys., Vol.169, pp.363-426, 2001.
11.25
1.5
1.75
Distance from origin[m]（ｙ＝０、ｚ＝０）6) 倉橋貴彦，斎藤浩一郎，近藤俊美，仮想領域法を用い2た 非圧縮性粘性流体の有限要素流体解析 － 定式化の方法の紹介と円柱周りの流れ解析 －， 長岡工業高等専

y=0，z=0 における x 軸上の温度分布（t=1.0sec.）

門学校紀要，Vol.50，pp.55-63，2014．

（非定常モデル：熱拡散率 k を変えた場合）

7) 日本数値流体力学会有限要素法研究委員会，有限要素法による流れのシミュレーション，シュプリンガー・フェアラーク東京，pp.58-59，1998．

５．おわりに8) 谷口健男，3 次元 FEM のための自動要素分割法，森北出版，pp.51-53，2006.11．

本論文では，三次元仮想領域FEMによる温度場解析に関する定式化を行い，数値実験による検討を行(2020. 8. 1 受付)った．定常モデルの解析を実施し，その結果を踏まえて，非定常モデルの解析を行った．非定常モデル
では，時間進行とともに定常解へ収束をし，定常モデルの解と一致することも確認した．また，熱拡散率を変えた場合においても，適切に解を得られていることを確認できた．将来的には，内部の埋め込み物体の領域の移動に伴った温度場の解析へ拡張を行う予定である．

謝辞：本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバ
SR16000を使用し計算を行ったものである．計算機の使用について，センター関係者の方へ謝意を表す．
参考文献
1) Ogasawara,K. and Kuroda,S.,Flow Analysis AroundOscillating Circular Cylinder Using Overset Grid Method,
14th Symposium on Computational Fluid Dynamics,
C01-2, pp.1-5, 2000.
2) Tang,H.S., Jones,S.C., Sotiropoulos,F., An overset-grid
method for 3D unsteady incompressible flows, J. Comput.
Phy., Vol.191, pp.567-600, 2003.

5