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

in 2 move

u  u の無次元速度の絶対値と無次元温度の分布をFig.4-2からFig.4-17に示す．円柱の壁面を白線で表し，図中のベクトルは速度を示す．

x y

u

move

u

 u

_in

u

_in

 u

_in

u

 u

_in

u

_in

 u

_in

L 2L 2L 2L L

L

0.5L 0.5L

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-2 u = u_in 2 _move^の結果（τ= 0）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-3 u = u_in 2 _move^の結果（τ= 1）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-4 u = u_in 2 _move^の結果（τ= 2）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-5 u = u_in 2 _move^の結果（τ= 4）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-6 u = u_in 2 _move^の結果（τ= 8）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-7 u = u_in 2 _moveの結果（τ= 1 2）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-8 u = u_in 2 _moveの結果（τ= 1 6）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-9 u = u_in 2 _moveの結果（τ= 2 4）

in 4 move

u  u の無次元の絶対速度と無次元温度の分布をFig.4-9からFig.4-15に示す．

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-10 u = u_in 4 _move^の結果（τ= 0）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-11 u = u_in 4 _move の結果（τ= 1）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-12 u = u_in 4 _moveの結果（τ= 2）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-13 u = u_in 4 _moveの結果（τ= 4）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-14 u = u_in 4 _move^の結果（τ= 8）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-15 u = u_in 4 _moveの結果（τ= 1 2）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-16 u = u_in 4 _moveの結果（τ= 1 6）

(a) 無次元速度の絶対値 (b) 無次元温度 Fig. 4-17 u = u_in 4 _moveの結果（τ= 2 4）

円柱が徐々に温められる様子を観察することができた． u_in2u_moveのとき円柱内の無次元最低温度は約0 . 8 2 9であり， 2 4において，u_in 4u_moveのとき円柱内の無次元最低温度は約0 . 9 2 4であった．u_in 4u_moveの方がu_in2u_moveに比べて入る熱量が多く，対流も強いため，乾燥炉内の温度が早く上がり，円柱の温度も上がりが早いことがわかる．また，円柱の動きにつられて全体の流れが右向きになった．流れが円柱を通りすぎる際に加速され，

in 2 move

u  u で無次元最高速度は約3 . 6 7，u_in4u_moveで無次元最高速度は約9 . 0 3となり，流入速度u_in^の2倍程度速くなった．そのため，ここで与えたR eに対して流れが不安定になった．u_in 2u_moveで圧力修正の反復回数は安定していたが，u_in 4u_moveでは周期境界付近で圧力修正の反復回数が収束しにくくなった．CPUにIntelのCore i7 3930K，GPUにAMDの Radeon R9 200 seriesを用いて計算した場合，u_in2u_moveは約1時間で計算することができ，

in 4 move

u  u は約50時間かかった．格子幅と時間刻みを更に細かくする必要があると考えられるが，GPUによる並列計算との相性がよく，実用的な手法であるといえる．

第 5 章

結言

本研究ではデカルト座標格子を用いて任意物体形状まわりの流れを精度良く計算するために，IB 法を応用し距離関数を用いて壁面温度境界条件を与える手法の導入とその手法の妥当性の検証を行った．

壁面で等温加熱条件を与えた計算の妥当性検証として，平板と格子のなす角度を変更し，

平板境界層流れの精度を確認した．平板と格子のなす角度が大きくなるにつれて精度が悪くなった．境界層内に格子が10格子程度用意すれば，速度と温度の分布は概ね一致することがわかった．また，壁面での速度と温度の勾配まで精度を安定させて解くためには境界層内に格子が20格子程度用意すればよいことがわかった．

物体内部の熱伝導を考慮し，対流の影響によって壁面での熱流束が異なる計算を行う手法を提案し，定性的に一致した結果を得た．

上で述べた手法を用いて，乾燥炉を模擬したモデル計算を行った．徐々に円柱が温められる様子を観測することができた．格子を十分細かくすることができれば，安定して計算することができ，実用的な手法であるといえる．

最後に今後の課題を述べる．

1. 物体内部の熱伝導を考慮した計算の定量的な評価

本研究では物体内部の熱伝導を考慮した計算の定性的な評価しか行っていないため，実験値や他の解析と定量的に比較する必要がある．

2. IB法の精度安定化

IB 法の特性として物体壁面と格子のなす角度によって精度が変わる^[9]．物体壁面と格子のなす角度が大きくなっても精度が悪くならない手法を導入する必要がある．

3. AMR (Adaptive Mesh Refinement)の導入

IB 法は物体壁面付近の格子を十分細かくする必要があり，計算領域全体を等間隔格子で計算することは効率が悪い．そのため，AMRを用いて物体壁面付近の格子を集中して細かくすることでメモリや計算時間を減らすことができ，効率よく計算することができる．

4. 乱流モデルの導入

自動車の塗装用乾燥炉では熱風の流入速度は1 0[m / s ]と大きく，自動車の大きさは4[m ] 程度あり，レールの移動速度は約0.1[m / s ]になるため，Re 2.0 10 ⁴となるため，LESなどの乱流モデルを導入する必要がある．

付録 A

支配方程式の無次元化

本研究における支配方程式は以下の3式である．

＜連続の式＞

  u (A.1)

＜Navier-Stokes方程式＞

 

¹ ²

u u u p u

t 



        

 (A.2)

＜エネルギー方程式＞

 

^u ²

   

     

 (A.3)

支配方程式の無次元化を示す．以下のように無次元変数を定義する．

X x

 x ,

U u

u ,

 t ,

P p

 p ,

0 c o l d

  



  (A.4)

式(A.4)を式(A.1)，(A.2)，(A.3)に代入し，整理するとそれぞれ次のようになる．ここで



^は

 

を無次元化したものであり，^{   }



^X^,^{ }^Y



^{と定義する．}

＜連続の式＞

 U = (A.5)

＜Navier-Stokes方程式＞

 

0 0

0 0 0 0 0

x U p

U U P U

u t u u x



 

       

 (A.6)

＜エネルギー方程式＞

 

0 0 0 0

x U

u t u x

    

 

     

 (A.7)

代表速度u₀は計算モデルによって一様流の流入速度や物体の移動速度とする．その他の代表値は式(A.8)のように定義した．

x0L, ₀_hot _cold, ⁰

0 0

x 1

u t  , ⁰

2 0

p 1

u ^ ^(A.8)

このとき未定参照量は次のように決まる．

0 0

t L

u , p₀u₀² _(A.9)

無次元数Reu L₀ ，Pr を用いて無次元支配方程式は次のようになる．

＜連続の式＞

 U = (A.10)

＜Navier-Stokes方程式＞

 

¹ ²

U U U P U

 Re

       

 (A.11)

＜エネルギー方程式＞





¹ ²

P r R e

  



     

 (A.12)

2.7.3節で述べた物体内の熱伝導を考慮した計算を行う際，エネルギー方程式は式(A.13)を

用いる．ここで，比熱Cは流体側では定圧比熱C_pを用いる．

 

^u ¹





t C

  



       

 (A.13)

熱伝導率kが一定のとき，kはナブラ演算子の外に出すことができる．このとき，k (C) より，式(A.13)は式(A.3)と等しくなる．上で述べた無次元化に合わせて無次元変数を式(A.14) のように定義する．

X x

 L,

U u

u ,

 L u , ^{co l d}

h o t co l d

  

 

 

 (A.14)

また，気相側の物性値を基準に取り，無次元の密度_C^，比熱C_C^{，熱伝導率}k_C^{を式(A.15)の} ように定義する．

C G

 

 , _C

C C

 C , _C

k k

k (A.15)

式(A.14)と式(A.15)を式(A.13)に代入し，整理すると次のようになる．

   

1 _G

C C G G

U k k

C C u L

  

  

      

 (A.16)

0 G

Reu L ^，Pr _G _Gを用いて無次元エネルギー方程式は式(A.17)のようになる．

 





C C

U k

C Pr Re

  

 

       

 (A.17)

付録 B

IP 点における値の補間

ある点 IPにおけるX方向速度U_IPの求め方について述べる．IP 点の座標



XIP^,YIP



が分かっているとする．始めに，その近傍の4 点のインデックスを求める．IP点の左下の近傍点のインデックスを



^{I J}^,



とすると，等間隔スタッガード格子において，



^{I J}^,



^{は式(B.1)，}

(B.2)で求めた



^I^{ }^, ^J



の商の整数部となる．

XIP

I  X

 (B.1)

IP 0.5

Y Y

J Y

 

   (B.2)

プログラミング上では式(B.3)，(B.4)のように整数型に変換することで，



^{I J}^,



^{を得ることが}

できる．

int X^IP

I X

 

    ^(B.3)

int Y^IP 0.5 Y

J Y

 

 

    ^(B.4)

例えば，



XIP^,YIP

 

 ^{1.08, 2.16}



，   X Y 0.1のとき，



^{I J}^,

 

^ ^{10, 22}



^となる．



^{I J}^,



^が求

まれば，



^I^^1,^J



^，



^{I J}^, ^¹



^，



^I^^1, ^J^¹



が残りの近傍点のインデックスとなる．Fig. B-1 のようにX_IP，Y_IPをとると，それぞれ式(B.5)，(B.6)のように求まる．

IP IP

X X I X

    (B.5)



^0.5



IP IP

Y Y J Y

     (B.6)

Fig. B- 1 IP点での補間

U

I_,

U

_I_₁_,_J

1 ,J

U

_I_₁_,_J_₁

XIP



IPY

U

) , (X_IP Y_IP

IP U

U

low

IP点の上方，下方の点におけるX方向速度をそれぞれU_up，U_lowとし，線形補間によって式 (B.7)，(B.8)のように求める．

1, 1 ( ) , 1

IP I J IP I J

X U X X U

U X

  

    

  (B.7)

1, ( ) ,

IP I J IP I J

low

X U X X U

U X

     

  (B.8)

Uup，U_lowを用いて，U_IPを式(B.9)のように線形補間して求める．

( )

IP up IP low

Y U Y Y U

U Y

    

  (B.9)

IP 点におけるY方向速度V_IPを求めるときはV の定義点の位置に合わせて式(B.3)，(B.4)は式 (B.10)，(B.11)のようになる．

int X^IP 0.5 X

I X

 

 

    ^(B.10)

int Y^IP

J Y

 

   ^(B.11)

また，式(B.5)，(B.6)も式(B.12)，(B.13)のようになる．



^0.5



IP IP

X X I X

     (B.12)

IP IP

Y Y J Y

    (B.13)

その後の操作はU_IPを求めるときと同様である．温度などのセル中心で定義される値に関しても定義点の位置に合わせてI，J，X_IP，Y_IPを求めることで，IP 点での値を補間することができる．

付録 C

境界条件の離散化

3.1 節の式(3.2)と式(3.3)で示した右壁と上面の境界条件の離散化について述べる．式(3.2) と式(3.3)を無次元化すると式(C.1)と式(C.2)のようになる．ここで，SとNは無次元の接線方向座標と法線方向座標であり，Ss L，N n Lと定義する．

右面： U 0 S

 

 ， 0



 

 (C.1)

上面： U 0 N

 

 ， 0



 

 (C.2)



^, ^,



f  U V  に対して f Sは式(C.3)のように表すことができる．

f f X f Y

S X S Y S

      

     (C.3)

cos

X S



   ，  Y S sin



より，式(C.3)は式(C.4)のようになる．

c o s s i n

f f f

S  X  Y

    

   (C.4)

f X

  と f Yを離散化して求めれば f Sの値を得ることができるが，右面の境界条件において， f Yを離散化することができない．そのため，Fig.C-1のY方向の勾配とX方向の勾配を用いて f Sを得る． f Yは式(C.5)のように表すことができる．

1 1

2 2

f f f

Y X Y

    

    (C.5)

また，式(C.4)と式(C.5)より，式(C.6)を得ることができる．



^{c o s} ^{s i n}



^{2 s i n}

f f f

S   X  Y

     

    (C.6)

f X

  と f Yを 1 次精度の片側差分法を用いて求めると f Sは式(C.7)のように表すことができる．ここで，  X Yとする．



^cos ^sin



^, ^1, ^{2 sin} ^, ^1, ¹

i j i j i j i j

f f f f

S   ^X^  ^ ^{ }X

   

   (C.7)

仮想格子内の値 f_{i j}_, に  f S 0となるように値を与えることで境界条件を与えることができる．つまり，式(C.7)に  f S 0を代入し，整理すると右面の境界条件は式(C.8)のようになる．

 

1, 1, 1

cos sin sin

cos

i j i j

i j

f f

f   



  

 

 (C.8)

36 また， f Nは式(C.9)のように表すことができる．

f f X f Y

N X N Y N

     

     (C.9)

sin

X N



    ，  Y N cos



より，式(C.9)は式(C.10)のようになる．

s i n c o s

f f f

N  X  Y

     

   (C.10)

f X

  と f Yを離散化して求めれば f Nの値を得ることができるが，上面の境界条件において， f X を離散化することができない．そのため，Fig.C-1のX方向の勾配とY方向の勾配を用いて f Nを得る． f Xは式(C.11)のように表すことができる．

1 1

2 2

f f f

X X Y

     

    (C.11)

また，式(C.10)と式(C.11)より，式(C.12)を得ることができる．



^{c o s} ^{s i n}



^{2 s i n}

f f f

N   Y  X

     

   (C.12)

f Y

  と f Xを 1 次精度の片側差分法を用いて求めると f Sは式(C.7)のように表すことができる．ここで，  X Yとする．



^cos ^sin



^, ^, ¹ ^{2 sin} ^, ^1, ¹

i j i j i j i j

f f f f

N   ^X ^  ^ X^{ }

   

   (C.13)

仮想格子内の値f_{i j}_, に  f N 0となるように値を与えることで境界条件を与えることができる．つまり，式(C.13)に  f N 0を代入し，整理すると上面の境界条件は式(C.14)のようになる．

 

, 1 1, 1

cos sin sin

cos

i j i j

i j

f f

f   



  

 

 (C.14)

Fig. C- 1 離散化の参照点（左：右面，右：上面）

45°



仮想格子

i, j i-1, j

i-1, j-1 X

Y Y’

135° 

i, j

i, j-1 i+1, j-1 仮想格子

X X’ Y

付録 D

Runge-Kutta 法を用いた Blasius 方程式の解法

3.1 節の比較に用いた値の求め方について述べる．ここで，一様流と平板はx軸に平行であり，平板の先端を原点とする．式(D.1)に示す平板境界層方程式はy uin

 

x の関数

 

F  を用いて式(D.2)のように変形することができる^[22]．^F

 

^ は式(D.3)のように定義され，

 

F   U u u となる．

2 2

u u u

u v

x y  y

    

   (D.1)

2FFF0 (D.2)

 

0 i n

F u d



   

 (D.3)

同様にして，式(D.4)に示す温度境界層方程式は無次元温度^{ }

 

を用いて，式(D.5)のように変形することができる^[24]．

u v 2

x y y

   

    

   (D.4)

2P r F0 (D.5)

ここではPrandtl数Pr0.7とする．また，境界条件は式(D.6)，式(D.7)の通りである．

 

⁰

 

⁰ ⁰

F F  ^，^{F  }

 

¹ (D.6)

 

⁰ ¹

  ，^{  }

 

⁰ (D.7)

4次精度のRunge-Kutta法は式(D.8)から式(D.12)のように求めることができる．Runge-Kutta

法を用いて式(D.2)と式(D.5)の解を求めるためには^F

 

⁰ ^と^

 

⁰ ^{の値が必要である．}^F

 

⁰

は Howarth によって求められた^F

 

⁰ ^^0.33206^[22]^{を用いる．}^

 

⁰ においては式(D.11)から

 

  ^が 0 になるように試行的に求め，^^

 

⁰ ^{ }^0.29268^{となった．このとき} を無限遠方とする．

    

1 2 3 4



1 2 2

F   F  6 k  k  k k (D.8)

    

1 2 3 4



1 2 2

F^  F^  6 l  l  l l (D.9)

    

1 2 3 4



1 2 2

F^  F^  6 m  m  m m (D.10)

    

1 2 3 4



1 2 2

6 c c c c

         (D.11)

ドキュメント内距離関数を用いた任意物体形状まわりの対流熱伝達計算 (ページ 32-58)