距 離 関 数 を 用 い た 任 意 物 体 形 状 ま わ り の 対 流 熱 伝 達 計 算

平成

30

年度（2018年度）学位論文（修士）

距 離 関 数 を 用 い た 任 意 物 体 形 状 ま わ り の

対 流 熱 伝 達 計 算

提出日：平成

31

年（2019年）

1

月

25

日

首都大学東京大学院

システムデザイン研究科 システムデザイン専攻 航空宇宙システム工学域 博士前期課程

学修番号

17891522

氏名 中谷 優浩

指導教員 田川 俊夫 准教授

i

目次

1.1

はじめに

1.2

格子の生成方法

1.3

塗装用乾燥炉の解析

1.4

本研究の目的

2.1

計算格子

2.2

支配方程式

2.3

無次元の支配方程式

2.4

離散化スキーム

2.5

距離関数

2.6 IB

法

2.6.1

滑りなし条件（No-Slip Wall）

2.6.2

滑り条件（Slip Wall）

2.7

壁面における温度境界条件

2.7.1

等温加熱条件

2.7.2

断熱条件

2.7.3

物体内の熱伝導計算

2.8

計算フロー

3.1

平板境界層流れ（等温加熱条件）

3.1.1

計算モデル

3.1.2

計算結果

3.2

対流熱伝達計算（物体内の熱伝導計算）

3.2.1

計算モデル

3.2.2

計算結果

第

1

章 緒言 ｐ.1

第２章 解析手法 ｐ.4

第

3

章 計算の妥当性検証 ｐ

.14

ii 4.1

計算モデル

4.2

計算結果

付録

A

支配方程式の無次元化

p.31

付録

B IP

点における値の補間

p.33

付録

C

境界条件の離散化

p.35

付録

D Runge-Kutta

法を用いた

Blasius

方程式の解法

p.37

付録

E  U

_S

 N

^と

N u

の求め方

p.40

付録

F

平板境界層計算の可視化結果

p.42

第４章 乾燥炉解析への適用 ｐ.23

第５章 結言 ｐ.30

付録

参考文献 ｐ.48

謝辞 ｐ.50

iii

記号表

A

：IB法の補間に用いる係数の分母

[  ] A

B

, A

_IP1

, A

_IP2 ：IB法の補間に用いる係数の分子

[  ]

C

C ^{：無次元比熱}

[  ]

C

G ：気相の比熱

[J/(kg  K)]

C

p ：定圧比熱

[J/(kg  K)]

C

S ^{：物体の比熱}

[J/(kg  K)]

C

^：比熱比

[  ]

C

：比熱

[J/(kg  K)]

d

：平行平板高さ

[m]

e

N ：無次元単位法線方向ベクトル

[  ] e

S ：無次元単位接線方向ベクトル

[  ] F

IB ：IB法による無次元外力項

[  ] H e

：Heaviside関数

[  ]

h

s ：局所熱伝達率

[W/(m

²

 K)]

k

C ^{：無次元熱伝導率}

[  ] k

G ：気相の熱伝導率

[W/(m  K)]

k

S ：物体の熱伝導率

[W/(m  K)]

k

^{：熱伝導率比}

[  ]

k

：熱伝導率

[W/(m  K)]

s l i p

L

：滑り面長さ

[m]

L

x ^{：計算領域の}

X

方向長さ

[m]

L

y ：計算領域の

Y

方向長さ

[m]

L

：代表長さ

[m]

N u

：局所ヌセルト数

[  ]

N

X ：

e

_Nの

X

方向成分

[  ]

N

Y ：

e

_Nの

Y

方向成分

[  ]

N

：無次元法線方向座標

[  ]

iv

n

：法線方向座標

[m]

P r

：Prandtl数

[  ]

P

：無次元圧力

[  ]

p

：圧力

[Pa]

Q

IB ：IB法による無次元熱量

[  ] R e

：Reynolds数

[  ] S

X ：

e

_Sの

X

方向成分

[  ] S

Y ：

e

_Sの

Y

方向成分

[  ] S

：無次元接線方向座標

[  ]

s

：接線方向座標

[m]

t

：時間

[s]

U

B ：物体壁面の無次元速度ベクトル

[－]

U

I B ：

IB

点における無次元速度ベクトル

[－]

U

IB ：

U

_{I B}の

X

方向成分

[－]

U

I P ：

IP

点における無次元速度ベクトル

[－]

U

IP ：

U

_{I P}の

X

方向成分

[－]

U

in ：無次元流入速度

[－]

N I B

U

：

U

_{I B}の法線方向成分

[－]

N I P

U

：

U

_{I P}の法線方向成分

[－]

U

N ：無次元法線方向速度

[  ]

S I B

U

：

U

_{I B}の接線方向成分

[－]

S I P

U

：

U

_{I P}の接線方向成分

[－]

U

S ：無次元接線方向速度

[  ] U

：無次元速度ベクトル

[  ] U

：無次元

X

方向速度

[  ]

u

in ：流入速度

[m/s]

u

move ^{：物体の移動速度}

[m/s]

u

：速度ベクトル

[m/s]

u

：

x

方向速度

[m/s]

v

V

IB ：

U

_{I B}の

Y

方向成分

[－]

V

IP ：

U

_{I P}の

Y

方向成分

[－]

V

：無次元

Y

方向速度

[  ]

v

：

y

方向速度

[m/s]

 X

：無次元

X

方向格子幅

[  ] X

：無次元座標ベクトル

[  ] X

：無次元

X

方向座標

[  ]

x

：座標ベクトル

[m]

x

：

x

方向座標

[m]

 Y

：無次元

Y

方向格子幅

[  ] Y

：無次元

Y

方向座標

[  ]

y

：

y

方向座標

[m]



：熱拡散率

[m

²

/s]



IP ：IP点までの距離

[  ]

1



IP ：IP1点までの距離

[  ]

2



IP ：IP2点までの距離

[  ]



IB ：物体近傍範囲

[  ]



：界面厚さ

[  ]



X ：距離関数の

X

方向勾配

[  ]



Y ：距離関数の

Y

方向勾配

[  ]



：距離関数

[  ]



：動粘度

[m

²

/s]



 ：物体壁面の無次元温度

[  ]



 ：IB点における無次元温度

[  ]



IP ：IP点における無次元温度

[  ]



：無次元温度

[  ]



hot ：最大温度

[K]



cold ：最小温度

[K]



：温度

[K]

vi



C ：無次元密度

[  ]



G ：気相の密度

[kg/m

²

]



S ：物体の密度

[kg/m

²

]



：密度比

[  ]



：密度

[kg/m

²

]





：無次元時間刻み

[  ]



：無次元時間

[  ]



：傾斜角度

[deg]

Notations

 

：

 

 





 

 













y x



^：

 

 





 

 













Y X



2 ^： 2

2 2 2

y

x 

 



 or

2 2 2 2

Y

X 

 





1

第

1

章

緒言

1.1

はじめに

CFD(Computational Fluid Dynamics)とはコンピュータによって流体に関する方程式を解き，

解析する分野である．計算機の発達とともに発展し，現在は航空機，船舶，自動車などの機 械設計や，血液の流れなどの解析といった生理学の解明，地球ダイナモや気象のシミュレー ションといった地球科学の解明や予想など多種多様な分野で用いられている^[1]．

CFD

は実験と比較すると設備投資のコストが低い，条件の変更が容易，実験での再現が 困難な現象の解析が可能，目視できない流れの可視化が可能などの利点がある．しかし，数 値解析によって求まる方程式の解は近似解であり，実際の現象を完全に再現することは困 難である．したがって，CFD を利用する場合には誤差がどの程度発生するのかを把握して 用いなければならない．

1.2

格子の生成方法

流体計算に用いられる格子は大きく分けて境界適合格子，非構造格子，デカルト座標格子 の

3

種類に分類される．それぞれの特徴を

Table1-1

に示す．デカルト座標格子は格子点あ たりの計算時間およびメモリが最小，格子生成が容易，空間高次精度が容易などの利点を持 つ．しかし，物体壁面が格子に沿っていない場合，物体形状の再現度を上げるために特別な 取り扱いが必要になる．また，十分な精度を確保するためには物体近傍で格子を細かくする 必要があり，必要なメモリが増加する．近年の計算機の発達によって，十分な格子数を確保 しつつ，実用的な計算時間で計算することができるため，デカルト座標格子の利用が見直さ れつつあり，複雑形状を対象にした流体解析への使用が期待される．

デカルト座標格子における物体境界の再現精度を上げる手法として

IB

法（Immersed

Boundary Method）

^[2]がある．IB法の優位性として以下の

4

点が挙げられる．

・格子生成にかかる手間が少ない

・物体形状が複雑化しても計算時間，計算負荷の変化が小さい

・移動物体境界問題への適応が比較的容易である

・高精度のスキームが利用できる

これらの利点から現在，盛んに研究が進められている．また，IB 法は格子と境界壁面のな す角度によって精度が変化するという問題がある^[3][4]．多くの

IB

法による解法が提案され ており，周囲の速度から補間して速度を与える手法[3][4][5][6][7]，物体壁面付近で壁面の位置を

2

考慮して離散する手法[8][9][10][11]などがある．また，物体を考慮して圧力修正の反復を行うこ とで精度を向上させる手法も開発されている^[12][13]．IB法を用いたモデルとして乱流解析な ど多くのモデルが検証されているが熱伝達を含むモデルは少なく，定量的に評価する必要 がある．移動物体まわりの熱伝達を計算することができれば，自動車の塗装用乾燥炉の解析 や食品加熱工程の解析といった加熱されながら物体が移動する生産ラインの最適化シミュ レーションに応用することができ有益であるといえる．

Table1- 1 格子の種類

^[14]

格子の種類 境界適合格子 非構造格子 デカルト座標格子

格子生成 ・簡単な形状は容易

・3 次元複雑形状には 単一格子の形成が困 難

・任意形状に対して自 動形成可

・解適合格子が容易

・格子作成の計算量が 大きい

・任意形状に対して自 動形成可

・解適合格子が容易

計算精度 及び 計算時間

・高効率，高精度な解 法

・物体表面近くの精度 が必要な問題に適し ている

・計算時間は境界適合 格子と遜色ないレベ ルに達した

・ メ モ リ 要 求 が 大 き く，空間高次精度化 も困難

・格子点当たりの計算 時 間 お よ び メ モ リ は最小

・物体境界の精度が特 別 な 扱 い な し で は 悪い

1.3

塗装用乾燥炉の解析

自動車の製造する際，車体を塗装液で満たされたプールに通し，車体表面をコーティング した後，左右から熱風を当てることで塗装液を乾燥させる処理がある．自動車製造において，

乾燥工程におけるエネルギー消費は全体の中でも大きな割合を占めており，その削減は重 要な課題である．

Fig1-1

に株式会社ディライトによって解析された例を示す^[15]．このとき車

3

1

http://ipconet.org/seminar/2014kyusyu/paper%205.pdf

体はレールに乗って乾燥炉内を等速で移動する．左右の側壁に熱風が出入りする流入口と 抜ける流出口があり，徐々に車体が熱される．現在，塗装用乾燥炉の解析には非構造格子を 用いることが多く，格子を形成するための計算負荷が大きく，メッシュを手直しする処理が 必要になることも多い．また，流体計算は計算に時間を要するため，熱風の流入口からの距 離や流入速度，流入角度などを用いて温度や熱伝達率をモデル化して解析されることもあ る^[16]．

IB

法を用いてデカルト座標格子で計算することで，格子生成にかかる手間を省くことが でき，プリ処理の負担を軽減させることができる．しかし，デカルト座標格子では非構造格 子を用いるより格子を細かくする必要がある．近年の計算機の発達によりデカルト座標格 子を用いて格子を十分細かくすることができ，かつ実用的な計算時間で解析することがで きる．今後，更なる計算機の発達によって移動物体の解析はデカルト座標格子が主流になる ことが考えられる．

Fig. 1- 1 塗装用乾燥炉の解析例

¹

1.4

本研究の目的

デカルト座標格子を用いて任意物体形状まわりの流れを精度良く計算するために，IB 法 を応用し，距離関数を用いて壁面温度境界条件を与える手法の導入とその手法の妥当性の 検証を本研究の目的とする．ここでは物体内部を等温とし，壁面で等温加熱条件を与えた計 算と物体内部の熱伝導を考慮し，対流の影響によって壁面で熱流束が異なる計算を行った．

4

章では応用例として乾燥炉を模擬した計算を行った．

4

第２章

解析手法

2.1

計算格子

本研究ではデカルト格子を用い，空間の離散化は

Fig.2-1

に示すような等間隔のスタッガ ード格子を用いた．圧力や温度などのスカラー量をセルの中心に定義し，ベクトル量である 速度をセルの境界に定義する^[17]．

2.2

支配方程式

支配方程式は以下に示す通りである．

＜連続の式＞

0

   u (2.1)

＜Navier-Stokes方程式＞

  ¹

²

u u u p u

t 



        

 (2.2)

＜エネルギー方程式＞

  ^u

²

t

   

     

 (2.3)

j

u

i_,

u

_i1,j

1, i j

u

_

2, i j

u

_

, 1

u

i j_ 1, 1

i j

u

_{ }

u

_{i j, 1}

1, 1 i j

u

_{ }

u

_{i j, 1}

,

p

i j , 1

p

i j_

, 1

p

i j_

1, i j

p

_

1, i j

p

_

,

v

i j

v

_i1,j

1, i j

v

_

, 1

v

i j_

, 1

v

i j_

v

_{i 1,j1}

1, 1 i j

v

_{ }

, 2

v

i j_

Fig. 2- 1 スタッガード格子の定義点

5

2.3

無次元の支配方程式

無次元変数と無次元数は以下のように定義する．

－無次元変数－

X x

 L ,

0

U u

 u ,

0

t

  L u ,

2 0

P p

 u

 ,

^{co l d}

h o t co l d

  

 

 

 (2.4)

－無次元数－

u L

0

R e ^  ^{, P r}



  (2.5)

無次元化した支配方程式は以下のようになる．導出は付録

A

を参考のこと．ここで，



は 無次元のナブラ演算子である．

＜連続の式＞

0

  U = (2.6)

＜Navier-Stokes方程式＞

  ¹

²

U U U P U

 Re

      

 (2.7)

＜エネルギー方程式＞

 ^U  ¹

²

P r R e

  



     

 (2.8)

2.4

離散化スキーム

本研究では時間項に式(2.9)に示すようなオイラー陽解法を用いた．

n n 1 n

U U U

 







  

   

  ^(2.9)

移流項に式(2.10)に示すような

UTOPIA

^[18]を用いた．UTOPIA では計算する点における速度 の正負によって風上を決定し，参照点に計算点と風上側に

2

点，風下側に

1

点の計

4

点を 用いて離散化する．計算領域の端と物体の近傍では

1

次風上差分法を用いた．

2 1 1 2

2 1 1 2

8 8

1 2

4 6 4

1 2

i i i i

i i

i i i i i i

U U U U

U U U

X X

U U U U U U

X

   

   

   

   

   

 

   

 

(2.10)

その他の離散化には

2

次精度中心差分法を用いた．

6

2.5

距離関数

距離関数とは壁面境界からの無次元符号付き距離であり，ここでは



とおく．本研究では 物体外部を正，物体内部を負とし，スカラー点とベクトル点のすべての位置で定義した．例 として，座標

  ^{1, 1}

^{を中心とする半径}

^0.5

の円柱を計算した．任意の点

 ^{X Y ,} 

^{における距離}

関数



は式(2.11)のように表すことができる．

 ^{X 1}  

²

^{Y 1} 

²

^{0 . 5}

      (2.11)

この円まわりの距離関数を可視化したものを

Fig.2-2

に示す．図中の黒線は

  0

の等高線で あり，円柱の壁面を示す．円柱の中心で最も値が小さく，計算領域の外側にいくに連れて値 が大きくなっていることがわかる．

Fig. 2- 2 円まわりの距離関数

本研究では距離関数を

IB

法の適用範囲の判定，

IB

法の補間に用いる係数の計算，物体と 流体の遷移領域の計算に用いる．また，距離関数



の

X

方向微分と

Y

方向微分をそれぞれ



X^，



_Yとすると壁面における無次元の単位接線方向ベクトル

e

_S，単位法線方向ベクトル

e

_N は式(2.12)のように表すことができる．これを用いて

IB

法における補間点の座標を計算と 速度の接線方向，法線方向への変換をすることができる．

2 2

2 2

Y

X Y

S

X

X Y

e



 



 

 

 

  

  

  

  

 

,

2 2

2 2

X

X Y

N

Y

X Y

e



 



 

 

 

  

  

 

  

 

(2.12)

7

2.6 IB

法

2.6.1 滑りなし条件（No-Slip Wall）

物体境界に隣接する格子の速度をどのように決定するかは様々な方法が考案されている．

本研究では計算の安定性が良く，補間が簡略である

Direct Forcing IB

法（直接強制法）^[3][4]

によって速度を埋め込む．

IB

法では式(2.13)のように

Navier-Stokes

方程式に外力項

F

_IBを加 える．

F

_IBは式(2.14)から式(2.16)のように物体遠方，物体近傍，物体内部によって異なる式 で表される．

U

_{I B}はまわりの速度から補間して求める速度であり，

U

_Bは物体の移動速度で ある．物体近傍の範囲を表す



_IB^は



_IB

 1.5  X

とすることで，物体壁面から少なくとも

1

格 子は物体近傍に含むことができる．

  ¹

² I B

U U U P U F

 R e

        

 (2.13)

物体遠方

 

IB

  

IB

0

F = (2.14)

物体近傍

 ⁰    

IB



  ¹

^{2 IB n}

IB

U U

F U U P U

Re 

       

 (2.15)

物体内部

 ^{  0} 

  ¹

^{2 B n}

IB

U U

F U U P U

Re 

       

 (2.16)

結局，物体遠方では

IB

法による操作を行わず，Navier-Stokes方程式の時間項にオイラー陽 解法を用いるため，物体近傍と物体内部では以下のように速度を与えることになる．

物体近傍

 ⁰    

IB



1 n

U

^

 U

IB

(2.17)

物体内部

 ^{  0} 

1 n

U

^

 U

B

(2.18)

U

Bは移動物体では物体の移動速度である．本研究では

U

_{I B}を以下のように与えた．

1 1 2 2

I P I P I P I P B B

I B

A U A U A U

U A

 

 (2.19)

     



^{1 2}1



² 1 ¹ 2



¹2



^{2 2 2 1 1}

, ,

IP IP IP IP IP IP IP IP IP IP

B IP IP IP IP

A A A

A

           

     

      

 

   

 (2.20)

8

式(2.19)は

Fig.2-3

に示すような

IB

法によって速度を与えている．これは物体壁面上の

B

点 の速度

U

_Bと物体壁面から



_IP1^離れた

IP1

点の速度

U

_IP1と



_IP2^離れた

IP2

点の速度

U

_IP2の計

3

点の速度を用いたラグランジュ補間によって物体近傍の

IB

点に速度を与えることを意味 する．

U

_IP1と

U

_IP2はそれぞれ

IP1

点と

IP2

点の近傍

4

点の速度を用いて線形補間し求め る．この補間方法に関しては付録

B

で詳しく述べる．また，



_IP1

 3.0  X

^，



_IP2

 4.5  X

^と した．

Fig. 2-3 物体近傍の補間（滑りなし条件）

2.6.2 滑り条件（Slip Wall）

滑り条件では物体遠方と物体内部は

2.6.1

節の滑りなし条件と同様にして与える^[3]．物体 近傍の速度

U

_{I B}は

Fig.2-4

のように壁面接線方向速度

U

_S^{と壁面法線方向速度}

U

_N^{に変換して} 条件を与える．

IB

点での壁面接線方向速度

U

_{S I B}と壁面法線方向速度

U

_{N I B}を求め，

U

_{S I B}と

N I B

U

を

X

方向速度

U

_IB^と

Y

方向速度

V

_IBに変換することで物体の影響を与える．壁面接線 方向速度は法線方向勾配 _S

0

U N

B



となるように与え，壁面法線方向速度は物体壁面で流 入出がないことが表されるように与える．

U

_{S I B}は式(2.21)のように物体壁面から



_IP^離れた

IP

点の壁面接線方向速度

U

_{S I P}をそのまま

U

_{S I B}に与える．これは

IB

点での法線方向勾配の

1

次精度前進差分が

0

であることに等しい．ここで，



_IP

 3.0  X

^とした．

S IB S IP

U  U (2.21)

N I B

U

は式(2.22)のように物体壁面上の

B

点の壁面接線方向速度

U

_{N B}と

IP

点の壁面接線方向

速度

U

_{N I P}の

2

点の速度を用いた線形補間で与える．

 

N IP IP N B

N IB

IP

U U

U   



 

 (2.22)

U

N Bは物体の移動速度を法線方向速度に変換して与える．本研究では静止物体でしか滑り 条件を用いないため，

U

_{N B}

 0

となる．

U

_{S I P}と

U

_{N I P}は

IP

点の速度

U

_{I P}を近傍の

4

点から補

U

IB

U

B 1

U

IP

1 IP

IB

B

IP2

U 2

IP

9

間して求め，壁面接線方向と壁面法線方向に変換し求める．単位接線方向ベクトル

 ^, 

S X Y

e  S S

，単位法線方向ベクトル

e

N

  N

X

^, N

Y



とすると，この変換は行列

A

を用いて 式(2.23)のように表される．また，IB点で求めた

U

_{S I B}と

U

_{N I B}は式(2.24)を用いて

U

_IB^と

V

_IB に変換する．

e

_S^と

e

_N^は

2.5

節で述べた通り距離関数



を用いて求める．

S I P I P I P X I P Y

N I P I P I P X I P Y

U U U S V S

U V U N V N

        

      

     

  A (2.23)

-1 S I B S I B Y N I B Y

I B

I B N IB S I B X N I B X

U N U S

U U

V U U N U S

  

 

       

       

  A     (2.24)

X Y

X Y

S S

N N

 

  

 

A _(2.25)

(a)

接線方向速度

(b)

法線方向速度

Fig. 2-4 物体近傍の補間（滑り条件）

IP

IB B

SIP

U

SIB

U

IP

IB B

U

N B

U

N IB

U

N IP

10

2.7

壁面における温度境界条件

2.7.1

等温加熱条件

速度境界条件と同様に

IB

法を用いて物体内部と壁面近傍の温度を与える．式(2.26)のよう にエネルギー方程式に熱量

Q

_IB^{を加える．}

Q

_IBは式(2.27)から式(2.29)のように物体遠方，物 体近傍，物体内部によって異なる式で表される．



_IBはまわりの温度から補間して求める温 度であり，



_Bは物体の壁面温度である．物体近傍の範囲を表す



_IB^{はここでも}



_IB

 1.5  X

とした．

 U  ¹

²

Q

I B

P r R e

  



      

 (2.26)

物体遠方

 

IB

  

IB

0

Q = (2.27)

物体近傍

 ⁰    

IB



  ¹

^{2 IB n}

Q

IB

U

PrRe

 

 



     

 (2.28)

物体内部

 ^{  0} 

  ¹

^{2 B n}

Q

IB

U

PrRe

 

 



     

 (2.29)

結局，物体遠方では

IB

法による操作を行わず，エネルギー方程式の時間項にオイラー陽解 法を用いるため，物体近傍と物体内部では以下のように温度を与えることになる．

物体近傍

 ⁰    

IB



1 n



^

 

IB

(2.30)

物体内部

 ^{  0} 

1 n



^

 

B

(2.31)



Bは物体の壁面温度であるため等温加熱条件では一定の値を持つ．本研究では



_IB^を以下 のように与えた．

1 1 2 2

I P I P I P I P B B

I B

A A A

A

  

  ^{ } (2.32)

     



^{1 2} 1



² 1 ¹ 2



¹2



^{2 2 2 1 1}

, ,

IP IP IP IP IP IP IP IP IP IP

B IP IP IP IP

A A A

A

           

     

      

 

   

 (2.33)

これは物体壁面上の

B

点の温度



_B^{と物体壁面から}



_IP1^離れた

IP1

点の温度



_IP1^と



_IP2^離れ た

IP2

点の温度



_IP2^の計

3

点の温度を用いたラグランジュ補間によって物体近傍の

IB

点に

11

温度を与えることを意味する．ここでも



_IP1

 3.0  X

^，



_IP2

 4.5  X

^とした．

2.7.2

断熱条件

断熱条件では物体遠方と物体内部は

2.7.1

節の等温加熱条件と同様にして与える．物体近 傍の温度



_IBは式(2.34)のように物体壁面から



_IP^離れた

IP

点の温度



_IP^{をそのまま}



_IB^に 与える．これは

IB

点での法線方向勾配の

1

次精度前進差分が

0

であることに等しい．ここ で，



_IP

 3.0  X

^とした．

IB IP

   (2.34)

2.7.3

物体内の熱伝導計算

エネルギー方程式のみ式(2.35)のように物体と流体の物性値の違いを考慮した式を用いる ことで，物体の熱伝導を計算する．比熱

C

は流体側では定圧比熱

C

_pを用いる．

  ^{u 1}  ^k 

t C

  



       

 (2.35)

無次元変数と無次元数は以下のように定義する．添え字

G

は気相側の物性値を表す．

－無次元変数－

X x

 L ,

0

U u

 u ,

0

t

  L u ,

2 0 G

P p

 u

 ,

^{co l d}

h o t co l d

  

 

 

 ,

C G

 

  ,

_C

G

C C

 C ,

_C

G

k k

 k

(2.36)

－無次元数－

0 G

R e u L

  ,

^G

G

P r 

  (2.37)

無次元化したエネルギー方程式は式(2.38)のようになる．

  ^{1 }

C

^

C C

U k

C Pr Re

  

 

       

 (2.38)

 k

c

 

  

は式(3.39)のように離散化した．

 

 

_,

,

1, , 1, , , 1, , 1,

, 1 , , 1 , , , 1 , , 1

2 2

2 2

C i j C C

i j

C i j C i j i j i j C i j C i j i j i j

C i j C i j i j i j C i j C i j i j i j

k k k

X X Y Y

k k k k

X X

X

k k k k

Y Y

Y

 



   

   

   

   

         

                    

   

  

 

   

  

 

(2.39)

物性値を物体と流体で区別した場合，格子形状に依存した結果になる．そこで本研究では物

12

体壁面に物体と流体の遷移領域を作り任意形状を再現した．遷移領域を持つ

Heaviside

関数

H e

^[19]を用いて各物性値を式(3.40)から式(3.42)のように計算した．添え字

S

は物体側の物性 値を表し，添え字

G

は気相側を表す．

 

G S G

He

       (2.40)

 

G S G

C  C  C  C He (2.41)

 

G S G

k  k  k  k He (2.42)

H e

は距離関数



を用いて式(2.43)のように表す． は遷移領域の厚さであり，

  1 . 5  X

と

した．

  X 1 . 0

のとき

H e

の形状は

Fig.2-5

のようになる．

 

 

 

1

1 1

1 s i n

2

0 H e

 

    

 

 

  

   

  

                

 



(2.43)

Fig. 2-5 遷移領域をもつ Heaviside

関数

式(2.40)，(2.41)，(2.42)の両辺をそれぞれ



_G，

C

_G^，

k

_Gで割ることで無次元変数



_C，

C

_C^，

k

Cは式(2.44)から式(2.46)のように表すことができる．ここで，



，

C

^，

k

^{はそれぞれ流体を} 基準とした物体と流体の密度比，比熱比，熱伝導率比である．

 

1 1

C G

 H e

 

     (2.44)

 

1 1

C G

C C C H e

 C    (2.45)

 

1 1

C G

k k k H e

 k    (2.46)

S G

 

 

， ^S

G

C C

 C

， ^S

G

k k

 k (2.47)

13

2.8

計算フロー

計算フローを

Fig.2-6

に示す．本研究では

GPU

による並列計算を行った．GPUを用いて 計算した処理は左に「GPU」と表記した．本研究では

Projection

法^[20]を用いたため，Navier-

Stokes

方程式によって速度の予測値を計算する際，圧力項を省いて計算し，圧力修正を行っ

た上で圧力項を足す．圧力修正の反復法は並列計算を行う際に隣接点の値のステップが混 同することを避けるために

Red & Black SOR

法^[21]を用いた．

Fig. 2-6 計算フロー

初 期 条 件 設 定

パ ラ メ ー タ 設 定

各種ファイル出力(可視化結果など)

速度の予測値を求める(圧力項以外)

IB法による速度補正

温 度 計 算

IB法による温度補正(等温加熱条件のときのみ)

圧力を計算する(Projection法)

速度の修正(圧力項計算)

物 体 移 動 ( 移 動 物 体 の と き の み )

時 間 進 行

各種ファイル出力(可視化結果など)

計 算 終 了 出力タイミング

時間終了 Yes

Yes

No

No GPU

GPU

GPU

GPU

GPU

GPU

GPU

距 離 関 数 計 算 距 離 関 数 計 算

GPU

GPU

14

第

3

章

計算の妥当性検証

3.1

平板境界層流れ（等温加熱条件）

3.1.1 計算モデル

IB

法の精度の確認ために壁面境界で等温加熱条件を与えた平板境界層流れを計算し，定 量的に評価した．格子に対する平板の角度

   0

，

1 5 

，

3 0 

，

4 5 

の

4

つ計算し，平板に 平行な一様流を与える．また，無次元格子幅を

  X 1.0 10 

^{2（代表長さに}

100

格子）,

5.0 10

3

X

^

  

^{（代表長さに}

200

格子）

, 2.5 10 

^{3（代表長さに}

400

格子）

, 1.25 10 

^3（代 表長さに

800

格子）の

4

つ計算した．Fig.3-1の(a)と(b)に示すように滑りなし条件と等温加 熱条件を与えた平板を長さ

3L

用意し，滑り条件と断熱条件を与えた平板をその前に長さ

L

用意した．

Navier-Stokes

方程式をもとに平板境界層流れを計算する際，滑りなし条件を与え た先端位置が特異点になり，そこから圧力勾配が発生する．この圧力勾配を解くために特異 点のまわりに空間を作る必要がある．そのため，流れに影響が出ないように滑り条件と断熱 条件を与えた平板を特異点の前に置いた．

   0

のとき法線方向の計算領域は平板から

L

だけ用意し，平板の高さ

d  0 . 1 L

とする．格子に平行でない平板の計算において，

Fig.3-1

の

(b)に示すように，少なくとも    0

の計算空間が含まれるように，

  15 

，

3 0 

，

4 5 

に

お け る 計 算 領 域 の

x

方 向 長 さ

L

_xと

y

方 向 長 さ

L

_yを そ れ ぞれ

 L

x

^, L

y

  ^{ 4.3 , 2.2} L L ^

，

 4. 1 , 3. 0 L L 

，

 ^{4.2 , 4.2 L L} 

^{とした．接線方向座標}

^s

^{と法線方向座標}

ⁿ

^{の原点は滑りなし条件}

と等温加熱条件を与えた平板の先端とする．境界条件は式(3.1)から式(3.4)の通りである．流 入面で平板に平行な一様流を与え，流出面では壁面接線方向勾配または壁面法線方向勾配 が

0

になるように設定した．圧力は全て反射条件とした．式(3.2)と式(3.3)の離散化と式展開 については付録

C

で述べる．

左面：

u  u

_in

cos 

^，

v  u

_in

sin 

^，

  

_cold ^{（物体外部）}

0

u 

，

  

_hot ^{（物体内部）}

^(3.1)

右面：

u 0

s

 



，

0

s



 

 (3.2)

上面：

u 0 n

 



，

0

n



 

 (3.3)

下面：

u  u

_in

cos 

^，

v  u

_in

sin 

^，

  

_cold ^{（物体外部）}

0

u 

，

  

_hot ^{（物体内部）}

^(3.4)

代表速度

u

₀

 u

_inとする．各パラメータを

Table 3-1

に示す．ここで，密度などの物性値は温

15

度依存しないと仮定し，浮力項は無視するものとする．

Table 3- 1 計算パラメータ

無次元格子幅

 X 1.0 10 

^2，

5.0 10 

^3，

2.5 10 

^3^，

1.25 10 

^3 無次元時間刻み

  1.0 10 

^4

Reynolds

数

1.0 10 

⁴

Prandtl

数

0 . 7

定量的な評価として，定常状態まで計算し，Fig.3-1に赤の点線で示した

s  L

の位置での無 次元接線方向速度

U

_Sと無次元法線方向速度

U

_N，無次元温度



の分布と壁面上の

 U

_S

 N

と局所ヌセルト数

N u

を用いた．本計算との比較に

Blasius

方程式を

4

次精度の

Runge-Kutta

法によって解いたものを用いた．この解法の詳細については付録

D

で述べる．また，

U

S

N

 

^と

N u

の求め方については付録

E

で述べる．

(a)

格子に平行な平板（

ψ = 0 °

）

(b) 格子に平行でない平板

Fig. 3- 1 計算モデル（平板境界層流れ）

Slip wall Thermal insulation

No slip wall Isothermal heating

L 2L

u

in

Observation line



in

d n s

1.1 L

L

Slipwall Thermalinsulation

Noslipwall Isothermalheating

L

2 L

Observation line

L

L

n s

 u

in



_in

L

x

y

L

16

3.1.2 計算結果

Fig.3-1

の赤点線部の

U

_S，

U

_N，



の分布を

Fig3-2

から

Fig.3-4

に示す．比較のために

Blasius

方程式を

4

次精度の

Runge-Kutta

法によって解いた値を黒の点線で示す．参考のために付録

F

で計算領域全体の可視化結果を示した．

(a) Δ = 1.0×10 X

^-2

(b) Δ = 5.0×10 X

^-3

(c) Δ = 2.5×10 X

^-3

(d) Δ = 1.25×10 X

^-3

Fig. 3-2 接線方向速度の分布

(a) Δ = 1.0×10 X

^-2

(b) Δ = 5.0×10 X

^-3

17

(c) Δ = 2.5×10 X

^-3

(d) Δ = 1.25×10 X

^-3

Fig. 3-3 法線方向速度の分布

(a) Δ = 1.0×10 X

^-2

(b) Δ = 5.0×10 X

^-3

(c) Δ = 2.5×10 X

^-3

(d) Δ = 1.25×10 X

^-3

Fig. 3-4 温度の分布

Fig.3-2

から

  X 1 . 0 1 0 

^2を除いて接線方向速度は概ね一致していることがわかる．

1.0 10

2

X

^

  

^では

   0

を除いて境界層が厚くなった．また，

  X 5.0 10 

^3の

  3 0 

，

4 5

において

N  7.0 10 

^{2付近で一様流より}

U

_Sが約

1%

大きくなった．

  X 2.5 10 

^3と

1.25 10

3

X

^

  

はよく一致した．Fig.3-3からすべての格子幅において法線方向速度は平板 角度によって異なる傾向を示した．特に，

  X 1.0 10 

^2の

  15 

，

3 0 

と

  X 5.0 10 

^3