第 25 回数値流体力学シンポジウム
E03-2

格子ボルツマン法を用いた電子デバイス直接冷却に関する液体流れの数値シミュレーション
Numerical Simulation of a Liquid Flow for Direct Cooling of an Electronic Device
Using a Lattice Boltzmann Method○ 高田尚樹, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1, E-mail: naoki-takada@aist.go.jp
ヨキネン アンテロ, VTT, FI-02150 フィンランド エスポー
コソネン トピ, VTT, FI-02150 フィンランド エスポー
松本純一, 産総研, 〒305-8564 茨城県つくば市並木 1-2-1
Naoki Takada, AIST, Namiki 1-2-1, Tsukuba, Ibaraki 305-8564, Japan
Antero Jokinen, VTT, Metallimiehenkuja 6, Espoo, P.O.Box 1000, FI-02150, Finland
Topi Kosonen, VTT, Metallimiehenkuja 6, Espoo, P.O.Box 1000, FI-02150, Finland
Junichi Matsumoto, AIST, Namiki 1-2-1, Tsukuba, Ibaraki 305-8564, Japan
A numerical simulation of an incompressible viscous fluid flow with heat transfer in a bending rectangular channel isconducted by using a lattice Boltzmann method, for evaluating basic performance characteristics of a small liquid
cooling unit directly attached on an electronics device. Heat conduction inside solid parts of the unit is also taken into
account. Major findings are as follows: (1) cooling rate is increased not only with the temperature difference between
the inflow fluid and the unit, but also with the flow rate; (2) for smaller thermal diffusivity, the cooling is slowed down,
but later, the rate becomes approximately constant for each thermal diffusivity; (3) the pressure at the inlet is increased
with the flow rate; (4) The unit is partly cooled down more rapidly around upstream side of the channel.

に様々な分野で熱・物質移動を伴う流体流れ問題(3)-(6)に適用されている．
2.2 格子ボルツマン法（LBM）
仮想的な
本研究の数値シミュレーションに用いる LBM(1), (2) は，粒子の粒子速度毎の存在確率（粒子数密度）fa の時間発展方程式を解くことで連続体近似の流体流れ場を予測する．
LBM の固体壁
境界はそこでの粒子の跳ね返り運動によって容易に設定できるた
め，LBM は複雑形状物体を含む流れ場の計算に適している(3)-(6)．
時刻 t，空間座標 x において，離散速度 ea (a: 速度の指標番号)を持つ粒子の速度分布関数 fa は次式に従って時間発展する．

１．はじめに近年，クラウドや携帯端末等を利用した情報・通信ネットワー
ク環境の高度化と各種サービスの拡大に伴って情報通信機器関連の消費電力量が増大している．集積度が進んだデバイス・機器・
データセンターの発熱は膨大であるため，それらの冷却に消費する電力を低減できるよう，室内空調を効率化する建築構造の改良
や，単相状態の液体や沸騰凝縮を伴う冷媒を利用した冷却システ
ムの開発と導入が国内外で進められている．本報では，上述のような Green IT の推進に寄与すべく，電子デバイス・機器が集積する各種装置や情報通信施設の冷却用消費電力の低減を目的とした，デバイスに直接載せる小型冷却ユニット内の液単相流れと固体および流体内部の熱移動に関する数値シミュレーションを述べる．

 fa
 e a   f a   a  x , t t

衝突による fa の時間変化a には BGK 近似 (1),(2)を適用した．

２．数値シミュレーション方法
2.1 計算対象
計算対象は，発熱する電子デバイス上部へ直接載せる小型冷却ユニット内部の，冷却用流路内を流れる非圧縮性単相 Newton 粘性流体および熱伝導を伴う固体部分である．流体は，ある一定流量で流入し，流路壁面から固体部分の熱を奪って外部の熱交換器へ流出，排熱することで電子デバイス冷却機能を果たす．その流体部分の圧力 p，流速 u および温度 T と固体部分の温度 T の分布の時間発展は次式に従って求める．
(1) u  0

u1
 u   u   p    2 u
t


(2)

T
 u  T   k  2Tt(3)

k  F , S (4) a ( x, t )  1 f ( x, t )  f aeq ( x, t ) f  a(5)

ここで平衡緩和時間は粘性などの輸送係数に関係する．単相流体の密度, 速度 u および圧力 p は次式で定義される．
(6)
   f a   f aeqaa u   a f a e a   a f aeq e a(7)p   cs2

(8)

上付添字 eq は平衡状態，a は全粒子に対する総和を意味する．cs
は音速である．
後述の計算では，一様な等方性空間格子上で 1 time step 間に運動粒子が隣接の格子点に到着する時間幅をt0 とし，式(4)を時間・
空間 2 次精度の Semi-Lagrange 形式(1),(2)に離散化した．
f a  x  e a t0 , t  t0   f a  x, t    a  x, t  t0 (9)ここで t は時刻， およびは流体の密度と動粘性係数，k は熱拡散係数であり，
下付添字kは流体( F )または固体( S )の値を表す．
本研究では上式を以下に述べる格子ボルツマン法（Lattice
Boltzmann Method, LBM）(1),(2)によって数値的に解く．LBM はすで

計算領域は，3 次元デカルト座標系(x, y, z)で単位長さの幅(x=y=z=1) の立方構造格子により一様に分割し，時刻を一定刻みで進行させる．粒子速度 ea については，静止状態(ea = 0 for a = 0)を含む 3 次元 15 速度モデル(a = 0～14)を採用した(1),(2)．今回使用
1

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2
2.4 非圧縮性流体流れのためのスキームの導入と改良
オリジナルの LBM の式(5)～(11)では，流れ場の代表速度 U が
音速 cs よりも十分小さい（低 Mach 数）条件下で非圧縮性流体流れが適切に再現される．このような LBM の弱い圧縮性がもたらす計算誤差は平衡分布関数faeqの式(12)の採用によって低減できる．
本研究では，式(12)を採用する代わりに，各時刻で非圧縮性の速度場u に関するソレノイダル条件を満たすように圧力p のPoisson
方程式を解くスキーム(9),(10),(11)を導入した．さらに，t0=1 以下の任
意の時間刻み幅t を設定できるよう時間進行方法を改良した．
改良スキームでは，各時刻，各格子点の粒子速度分布関数 fa お
よび ga,，流れ場の変数 u, p および T を以下の手順(1)～(6)に従って計算する．
(1) 現在時刻 tn の各格子点で衝突および並進演算を行う．

した ea=(ea,x, ea,y, ea,z)T は次式で与えられる．

e0 , e1 , e2 , e3 , e4 , e5 , e6 , e7 , e8 , e9 , e10 , e11 , e12 , e13 , e14  0 1 1 0 0 0 0 c  0 0 0 1 1 0 0
 0 0 0 0 0 1 1

1 1
1 1

1 1 1 1 1
1 1 1 1 1 
1 1 1 1 1 1 1 1 
1
1

(10)

LBM の平衡分布関数 faeq は，fa の時間発展式(1)から流体運動の支配方程式を漸近理論展開に従って導出できるように規定され，気体分子に関する Maxwell 分布から低 Mach 数の条件で流速 u に関する Taylor 展開することによって導出される(1),(2)．平衡分布関数
は以下で表される．

 e  u  e  u 2 u  u 
f aeq  wa  1  a 2  a 4  2 cs
2cs
2cs 


f a '  f a ( x , t n  t 0 )

(11)



式中の重み係数 wa には，静止粒子(a=0)に対して 2/9，運動粒子に対して 1/9 (a=1～6)，1/72 (a=7～14)が適用される．式(11)は圧縮性
=0 =constant の非圧縮
流れの支配方程式をもたらすのに対して，性流れの式(2)は次の平衡分布の適用によって導出される(7),(8)．

f

eq
a


 e  u  e  u 2 u  u  
 wa    0  a 2  a 4  2  
2 cs
2cs  

 cs





  cs2  f 

fTu   a g aeq e a(16)

 g a  x, t  

t0 g ,k g a  x, t   g aeq  x, t  a

t 0 2 

k  F, S 

(25)

0
t

(26)

  u ''

本研究では，上式を 2 次精度中心差分で離散化し，SOR 法を使用した収束演算により p n+1 を決定した．
(4) 時刻t n+1 で連続の式(1)が満たされるように圧力勾配でu''を修正して流速 u n+1 を求める．

(17)

u n 1  u ''

t

0

(27)

p n 1

右辺の圧力勾配も 2 次精度中心差分で離散化する．
(5) 分布関数 fa を次式によって更新する．

f a (x, t n 1 )  f a (x, t n  t )(18)

 f a '(19)


wa e a  t1 a
  p n 1 
 u '' un 
cs2  0a(28)(6) 温度 T とその分布関数 ga を以下のように更新する．

g a '  g a ( x, t n   t 0 )は，次式で定義されるTの拡散係数k の調整パラメータである．




(24)

t c  tt 0 x

 2 p n 1 2
ea  u  ea  u  u  u  2
g  x, t   Twa 3  2 cs
2cs4
2cs  k  c 2   g , k a(3) 新しい時刻 t n+1 =t n+t の圧力 p n+1 を以下の Poisson 方程式を解いて求める(9),(10),(11)．

eqau u 
g 0eq  x, t   T 1  3 1  w0   w0 2 
2cs 


(22)

a は任意の時間刻み幅t と格子単位時間幅t0 との比であり，格子定数 c に関する Courant 数に相当する．

下付添字 k は流体(F)または固体(S)の識別子である．本研究では関数 gaeq を次式で与えている．

 a  0au ''  a u ' 1  a  u n

T の時間発展式(3)も，式(5)および(9)と同じ形式で解く．

g a  x  e a t0 , t  t0 

 e  u  ea  u 2 u  u 
 wa 1  a 2 
 2
cs
2cs4
2 cs 
a

2.3 温度に関する計算
温度 T は，分布関数 ga を用いて式(15)で定義される．その平衡分布関数 gaeq は T と u との関係式(16)を満たす．
(15)

eq
a



(14)

T   a g a   a g aeq

 f (x  e a t0 , t n )  f aeq (x  e a t0 , t n ) f  a(2) 次の時刻 t n+t の仮の流速 u''を，現時刻 t n の un と並進後の関
数 fa'から次式で求める．これは u''を前後の時刻 t n と t n+t0の流速から線形内挿によって予測することを意味する．
(23)
u' 
f 'e

(12)

(13)

t 0 

2 

t 0

ここでの平衡分布関数は以下のように定義される．

上記の空間・時間・粒子速度の離散化により模擬される流体中の圧力p，
音速cs と動粘性係数 は以下のように導かれる (c =x/t0)．

cs  c / 3(21) f a ( x  e a t 0 , t n )(29)

 g a ( x  e a t 0 , t n )

(20)


2

t 0

 g ( x  e a t0 , t n )  g aeq ( x  e a t0 , t n ) g  aCopyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2(ii) 連続流出境界：

g 0 (x, t n  t )

(30)



 g 0 ' 1  a  T ' T n  1  3  wa a0p  p0  const. ， n u  0 ， n T  0

３．Benchmark Test Simulation
本研究では，温度 T に関する LBM の計算精度を簡潔に評価するため，
分布関数 ga の時間発展方程式(17)および平衡分布関数(18)，(19)を用いてFig. 1に示すような2種類の一次元温度拡散問題を解析した．初期条件では，赤色部分に T=T0=0.01，他の部分には T=0を設定した．3 次元計算領域の x 軸方向の両端には，一様な静止条件(a) u=0 では T=0 の Dirichlet 境界条件を適用し，一様に x 方向のみ定常移流する条件(b) u = u0 では周期条件を適用した．他の y, z
両方向には条件(a)，(b)ともに周期境界を設定した．なお，本解析
では，fa の時間発展を計算せずに各格子点で u を所定の値に規定
し，時間進行はオリジナルの時間刻み幅t0(=1)で行った．

g a ( x , t n  t )
e u
 g a ' wa 1  a  T ' T n   3  a 2
cs


n 1





 a  0

(31)

時刻 t n，仮の時刻 t n +t0 および次の時刻 t n+1 = t n +t における
温度 T は各々以下のように定義される．
(32)
T n ( x )  g ( x, t n )



a

a

T '(x)   g a '(x)(33)

T n 1 (x)   g a (x, t n  t )(34)

Periodic boundaries
a=20x

a

z

y

T=0

a

z

y

u0

1

DimensionlessTemperature T/T0

0.8

(37)

k  F, S 

0.6

T0=0.01
Initial value(at t = 0)

0.4

FDM
(2nd order)

Results by LBM
t*=16.7
t*=33.3
t*=66.7
t*=133.3
t*=266.7

g=1,  =1/3
t*=t/x2

0.2

(39)

0.4

x/Lx

0.6

0.8

1

(a) The numerical results under the stationary condition
(40)

1

(41)

流入・流出境界に関しては，計算領域外部から境界上の格子点へ進入する粒子の速度分布関数は，平衡分布関数 faeq および gaeq
によって与えた(12)．計算領域外部の境界に隣接する格子点での圧，u および温度 T は以下の条件から設定した．な
力 p（または）お，勾配は 1 次精度で差分近似した．
(i) 一様流入境界：n p  0 ， u  uin ，T  Tin

T   2Tt
=1.6710-1

0
0

DimensionlessTemperature T/T0

qW t k ck xb

Lx=256x

0.2

ここでx は格子点間距離，T-1 および T0 は壁面隣接の固体および流体格子点での温度である．熱伝達による流体および固体格子点での温度変化TF とTS を次式で求め，T-1 および T0 に加える．

Tk 

b=5x

T=0

(b) Advection condition
Fig. 1 One-dimensional diffusion problems.

(38)

 F T0   S T1
F  S

Periodic boundaries

x xC

ck は比熱，k は固体（k=S）または液体密度（F =0）である．本
研究では式(37)を一次精度で離散化して壁面温度 TW を決定する．

TW 

Lx=100x

Periodic boundary
a=35xT0=0.01OS およびF は固体および流体の熱伝導係数である．

T T
T T
 S W 1   F 0 W
x / 2
x / 2b T=0(a) Stationary condition(36)
u ''  0
ここで n は境界面の単位法線方向ベクトルを表す．
流体-固体間で熱伝達がある場合は，壁面上で熱流束 qW が連続であると仮定した．

 k   k ck k

b=5x

T=0

x

(35)

qW   S  n T  solid   F  n T  fluid

T0=0.01

T=0O

2.5 境界条件
本研究では，固体物体表面や計算領域境界面は格子点間の中間に位置すると想定する．流体領域に接しかつ固体領域内に位置する格子点では，分布関数 fa および ga を進入してきた方向へ速度-ea
で跳ね返す Bounce-back 条件を適用し，滑りなし静止断熱壁境界を再現する．改良スキームを採用して圧力 Poisson 方程式を解く
際は，壁面上で法線方向圧力勾配0 と仮流速u''0 の条件を課した．

n p  0

(43)

Results by LBM
g=1,  =1/30

t*=16.7
t*=33.3

0.8

 =1.6710-2

t*=66.7

T0=0.01
C=u0t/x=0.02
Pe=u0x/=1.2

t*=133.3

0.6: Results
for u0=0

0.4
t*=266.7

0.2
0
-0.2

-0.1

0

(x-xC-tu0)/Lx

0.1

0.2

(b) The numerical results under the advection condition
Fig. 2 Time evolutions of temperature distribution in the
one-dimensional diffusion benchmark problems.

(42)
3

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2
時間発展する T の分布の数値解析結果を Fig. 2 に示す．図中の○等のシンボルが LBM による結果を示し，黒線は式(2)を空間・時間とも 2 次精度の有限差分近似で陽的に直接解いて得られた結果である．両数値解の比較から，静止条件(a)，一様定常流条件(b)でともにLBMが温度の拡散を適切に予測することが確認できる．

Leak-proof rings(rubber)

Gap
space

４．冷却ユニット内熱移動・流体流れの数値シミュレーション
4.1 計算対象
電子デバイス上に直接載せる小型液体冷却ユニットを模擬した
計算対象をFig. 3 およびFig. 4 に示す．
3 次元デカルト座標系の下，計算領域を x, y および z 軸方向に各々232，232，71 個の格子点（幅x=y=z=1 の単位立方セル）で等分割し，冷却ユニットの実際の大きさ 46mm×46mm×13.8mm に対して空間解像度は x, y, z 各方
向 0.2mm とした．領域内部には，Fig. 3 に示されるような一続きの矩形断面流路が配置されている．冷却液に相当する非圧縮性流体は，一辺 d = 20x の正方形断面の開口境界からある一定の温度
Tin と速度 Uin で一様流入し，断面積の異なる流路内を流れながら，最終的には正方断面積 d 2 の流出境界から流出して一定圧力
p0=ρ0cs2 の下に解放されるとした． 重力の影響は無視した．
計算領域周囲には，厚さ 1 格子点分（1 セル幅）の断熱固体壁
条件を適用した．そのすぐ内側下部には面積 224x×224y，厚さ
1 格子点分の発熱固体層を，さらにその上部にはユニット部品間
ギャップ 1 か所および液漏れ防止用リング状ゴム材 2 点を模擬した断熱固体壁 3 点を配置した．上記以外の領域は熱伝導性固体が詰まっている．
実際の冷却ユニットの周囲には空気が存在するが，本研究では液体流れによる熱移動と電子デバイス冷却に注目するため，デバイスと外部空気との熱伝達は無視して断熱条件を周囲境界に課した．また，冷却ユニット下面は発熱する電子デバイスに密着していることから，発熱固体層はデバイスと同じ単位時間
あたり発熱量を持つと想定した．
本研究で想定した固体および液体の物性を Table 1 に示す．固体はアルミニウム合金（EN AW-6060），液体は密度，粘性，比熱についてはエチレン-グリコール(EG)重量比 40%水溶液（1atm, 10C）のものに相当する．
熱伝導および温度拡散については EG40%w 水溶液よりも約 17.7～70.5 倍大きい値に設定した．EG 水溶液の粘性は温度によって異なるため，本研究でも想定する流入流体温度に応じて動粘性係数の値を変更した．ただし，本研究は当該冷却ユニットの基礎的性能評価の第一歩であること，また，熱伝導性に関しては温度依存性が小さいことから，今回のシミュレーショ
ンではこれら物性は流れ場内の温度変化に依らず各流入条件下で常に一定と設定した．

0.85d
0.75d

7.1d
3d

Solid



Inlet

2.35d

d 1.55dd=20xx=y=z=10.75d

1.058×10(EG40%wt,
10C)

3

2.7×103[Kg/m ]

2.552

3Specific heatc
[J/(Kg K)]

3.621×10(EG40%wt,
10C)

2.458×10

2.0×107.888~31.34-1(cf. 4.448×10
for EG40%wt,
10C)

6.382~25.35
2(cf. 4.496×10
for EG40%wt)

8.323×10

2.059×10-6~8.181×10-7(cf. 1.161×10
for EG40%wt,
10C)

2

8.90×10

Thermalconductivity2

[W/(m K)]

-6

Thermaldiffusivity =  / c
2
[m /s]-5-1

10.17~40.43
2(cf. 7.169×10
for EG40%wt)

本計算対象の流れ場を代表する無次元数として，以下のReynolds 数 Re，Prandtl 数 Pr および Peclet 数 Pe を選択する．

d U in ,Pr 

 , Pe  Re  Pr
F

(44)

シミュレーションは，Re = 51.2～172.1，Pr = 1.0～4.0，Pe = 51.2
～319.3 の範囲の条件下で実施した．これは，実際の流入境界の高さと幅 d = 4mm と Table 1 の物性に対して体積流量 Qin = 0.1～0.3
LBM では，l/min., 流入速度 Uin = 0.1042～0.3125 m/s に相当する．
流体密度 ρ0 =1，代表長さ d = 20, 一様流入速度 Uin=5×10-3，時間刻
み幅t=0.5 を全条件で共通に使用した．また，分布関数 ga の計算
では，以下のように定義される無次元温度 T *を使用した．

0.85d
3.45d
0.6d

3d

d

Ratio
(Solid/Fluid)

Fluid
3

Density

0.75d

1.25d

OutletTable 1 Physical properties of solid and fluid simulated in this study.

7.8d

2.25d
0.35d

d

Fig. 4 Thermal conditions of the solid part.

11.5d

ChannelHeat-insulating
inside parts

Bottom heated
layer

Re 

Computational domain:
232x×232y×71zdHeat-insulating
outside wall

7.1d

T* 

11.5d

1.05d
: Main stream
direction

Fig. 3 The computational domain for simulation of heat transfer andfluid flow in a liquid cooling unit added on electronic device.
4

T  Tin
TS ,0  Tin(45)

ここで TS, 0 は初期固体温度，Tin は流入流体温度を表す．発熱固体
層では，デバイスの発熱による温度上昇T=(dT/dt)t を実際の冷却ユニットの熱設計値dT/dt =9.7C/sに基づき各時刻でTに加えた．
初期固体温度 TS, 0 は内部で一様とする一方，流路内の流体初期
温度 TF,0 は，以下の 2 種類の初期条件（I.C.）いずれかに設定した．
Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2
Fig.6 は，初期条件 Type(A)，流入流体温度 Tin=-10C で Pr と流量 Qin の異なる条件における，固体部分の平均温度 TS の時間変化を示す．
流量が多い場合固体はより速く冷却されることが分かる．
また，より大きな Pr=2 の場合，TS の減少は Pr=1 の場合より緩やかになる．これは，流体の温度拡散係数F がより小さいために固体境界付近から流体内部へ広がる熱量が減少するためである．

I.C. Type (A): TF , 0 = TS, 0
I.C. Type (B): TF , 0 = Tin (<TS, 0)

Solid temperature TS (C)

4.2 結果と考察
本研究で得られた流路水平断面内流速分布および流路壁面上の圧力分布を Fig. 5 に一例として示す．流れは Re=153.6 で定常状態に達している．図内のカラースケールには LBM の値と Table 1 の物性を持つ流体に対する実数値を示している．図(a)からは流体が流路に沿って蛇行しながら狭い流路断面で加速していることが確認される一方，図(b)では流入側から流出側に向かって流路に沿って圧力が次第に減少していることが分かる．

30

Tin = -10 deg.C

Pr = 2.029

Qin = 0.15 l/min.

28
27
Qin = 0.30 l/min.

260

Pr = 1.0

0.2

0.4

0.6

0.8

1

Time (s)Fig. 6 The cooling rates dependent on Qin and Pr, for Tin= -10C andI.C. Type (A).

u [m/s]2.25

Solid temperature TS (C)

0.0

u (LBM) [-](a) Flow velocity

30

Pr = 2.0
Qin = 0.15 l/min.

29
Tin = 0 deg.C

28
Tin = -10 deg.C

27
26
0

0.2

0.4

0.6

0.8

Time (s)1(a) Pr = 2.0 and Qin = 0.15 l/min.

p [MPa] 0.05

Solid temperature TS (C)

0.0

p (LBM) [-](b) Pressure on channel wallFig. 5 The numerical results for Re =153.6 (Qin= 0.3l/min).
以下に示すシミュレーション結果は，低 Re 数条件でほぼ定常状態に達した流速分布を初期条件に使用して得られたものである．
また，いずれの流動条件においても，固体部分の初期温度は実在のアルミニウム合金で 30C と想定した．
Figs. 6～10 は，流量 Qin，流入流体温度 Tin，Prandtl 数 Pr および
流体初期温度に関する各条件で得られた，熱伝導性固体部分の平均温度の時間変化を示す．図中の実在系の時刻(s)は，次の無次元
時刻 t*を通して LBM の値から求めた．

t* 

U in
t
d

30

Pr = 4
Qin = 0.10 l/min.

29

Tin = 0 deg.C

28
Tin = -10 deg.C

27
26
0

0.5

1

1.5

Time (s)(b) Pr = 4.0 & Qin = 0.10 l/min.
Fig. 7 The effect of inflow temperature Tin on cooling rate.

(46)

Fig. 7 (a)および(b)は各々，Pr と流量 Qin は同じで流入温度 Tin が5

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2
度 TS の時間変化を Fig. 9 (a)および(b)に示す．流入温度および流量は同じである． 両図を比較すると，流動開始 0.2 秒後までは Type(B)の方が TS はより急激に減少するが，その後 TS は両初期条件でほぼ同じ一定割合で減少することが確認できる．また，Pr 数に対する冷却速度の変化の傾向も両条件で同様である．
Fig. 10(a)および(b)は，初期条件 Type (B)，Pr=2 に対する，流入温度Tin と流量Qin を各々変更した場合の固体温度TS の時間履歴を表す． 初期条件 Type (A)に対する結果（Fig. 7）同様，Tin の低下と Qin の増加はより速い固体の冷却に繋がることがわかる．

30

Pr = 4.0
Tin = -10 deg.C

30

Solid temperature TS (C)

Solid temperature TS (C)異なる場合の固体温度 TS の時間変化を表す．両図から，より低温の流体の流入により固体はより速く冷却されることが分かる．
本シミュレーションで最も大きなPr 数条件で流量Qin を変化させた場合の固体温度の時間変化を Fig. 8 に示す．上記の Fig. 6 で示された結果と同様，
流量が減少すると冷却速度も減少している．
このことから，EG 水溶液のようにより大きな Pr の条件でも冷却速度には流量の変化に対する同様の傾向が現れると考えられる．

29
Qin = 0.10 l/min.

28
27

Qin = 0.15 l/min.

26
0

0.5

1

1.5

Time (s)

30

28

26
0.2

0.4

0.6

Time (s)0.8

Tin (deg.C):-10 (Re = 76.8)
24
0 (Re = 119.7)
5 (Re = 144.8)
10 (Re = 172.1)

22
0

1

(a) Initial condition Type (A)Solid temperature TS (C)

Pr = 1.0
Pr = 2.0
Pr = 4.0

0.4

0.6

Time (s)0.81

Tin = 0 deg.C
Pr = 2.0
Qin (l/min.):
0.15 l/min.
0.20 l/min.

24
22
0

0.2

0.4

0.6

Time (s)0.81

上記の各条件で観察された最大圧力増加Pmax を Table 2 に示す．
同じ流量では Re 数が低い
（低
Pmax は流量の増加とともに増加し，温で粘性が大きい）場合により大きくなることが分かる．
Table 2 Maximum pressure increase Pmax (measured at inlet corner).
Qin [l/min.] Tin [deg.C] Uin [m/s]

0.2

0.8

(b) Tin=0C and Qin= 0.15 and 0.2 l/min.
Fig. 10 Time history of solid temperature for Pr = 2.

24
22
0

0.6

26

Tin = -10 deg.C
Qin = 0.15 l/min.

26

0.4

28

30
28

0.2

30

27

0

26

(a) Qin= 0.15 l/min and Tin=-10C to +10C.

Tin = -10 deg.C
Qin = 0.15 l/min.
Pr = 1.0
Pr = 2.0
Pr = 4.029

28Time (s)Solid temperature TS (C)

Solid temperature TS (C)

Fig. 8 The cooling rate dependent on Qin for Pr = 4.0 and Tin= -10C.

Pr = 2.0
Qin = 0.15 l/min.

0.10
0.10
0.15
0.15
0.15
0.15
0.20
0.30

1

(b) Initial condition Type (B)Fig. 9 The cooling rate dependent on Pr, for Qin=0.15 l/min., Tin= -10C.

-10
0
-10
0
5
10
0
-10

0.10420.10420.15630.15630.15630.15630.20830.3125Re

Pmax* [-] Pmax [MPa]51.19479.82576.792119.738
144.759
172.111159.625
153.583

6.054E+02
5.709E+02
5.722E+02
5.623E+02
5.614E+02
5.251E+02
5.467E+02
5.459E+02

6.951E-03
6.554E-03
1.478E-02
1.452E-02
1.450E-02
1.356E-02
2.511E-02
5.641E-02

最後に，計算領域内部の水平断面内温度分布を Fig. 11(a)～(f)に
示す．
本結果は，
初期条件 Type (A)，
物性および流入条件 Pr=2, Qin=

流体温度の初期条件 Type (A)および Type (B)での平均固体温6Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
E03-2能を有することを確認した．また，液体による電子デバイスの直接冷却の効率化には，冷却液の温度と流量設定が重要であることも分かった．
当該冷却ユニットを含む液体冷却システムにより電子デバイ
ス・システム機器・設備を効率的に冷却するためには，目標とす
るデバイス温度，外部に接続するポンプの駆動エネルギー，熱交換器での排熱量，冷却ユニットと接続配管を含むシステム全体の圧力損失を適切に設定し，これら流体機器を設計・統合する必要がある．また，大型機器やそれらが集積する設備の低消費電力での冷却の実現には，エアコン装置等による全体的な空冷と局所的な液体冷却のベストミックスも不可欠である．

0.15 l/min., Tin= -10C の下，流動開始後の時刻 0.848s におけるものである．固体温度は上流側でより大きく減少し，その低温領域は上流側に偏って流路周辺に広がっている．また，流路内の温度の空間変化から，流体が流路壁面を通して固体から熱を奪って輸送していることも確認される．このことから，本計算対象では，固体は流路に沿って上流側から優先的に冷却され，流体-固体間の熱伝達で温度の低下した流路周辺へ向かってより内側の固体部分の熱が熱伝導により移動して冷却が進展すると考えられる．

T-Tin
(C)
40

謝 辞
上記の研究成果は，VTT（Technical Research Centre of Finland, フィンランド国立技術研究センター），Tekes（The Finnish Funding
Agency for Technology and Innovation, フィンランド技術庁），フィ
ンランド国内民間企業数社，および産総研の共同研究プロジェク
トにおいて得られたものである．関係者各位への謝意をここに表します．

0(a) Height z=0.1mm(d) z=3.9mm(b) z=1.3mm(e) z=8.9mm(c) z=2.5mm(f) z=13.7mm

参考文献(1) Chen, S. and Doolen, G. D., "Lattice Boltzmann Method for Fluid
Flows," Annu. Rev. Fluid Mech., 30 (1998), pp.329-364.
(2) Succi, S., "The Lattice Boltzmann Equation for Fluid Dynamics andBeyond," (2001), Oxford at the Clarendon Press, UK.
(3) Yoshino, M. and Inamuro, T., "Lattice Boltzmann simulations offlow and heat/mass problems in a three-dimensional porous
structure," Int. J. Numer. Meth. Fluids, 43 (2003), pp. 183-198.
(4) Yoshino, M., Murayama, T., Matsuzaki, A. and Hitomi, T., "Masstransfer analysis of calcium in concrete using the lattice kinetic
scheme for a binary miscible fluid mixture," Transactions of theJapan Society of Mechanical Engineers, Series B, 73 (2007), pp.
973-980.
(5) Yamamoto, K., Yamauchi, K., Takada, N., Misawa, M., Furutani, H.,
Shinozaki, O., "Lattice Boltzmann simulation on continuously
regenerating diesel filter," Philos. Trans. R. Soc. A-Math. Phys. Eng.
Sci., 369 (2011), p. 2584-2591.
(6) Takada, N., Hayashi, H., Seta, T., Yamamoto, K., Yoshino, M. andMatsukuma, Y., "Numerical Analysis of Complex Fluid Flows by
the Lattice Boltzmann Method," Journal of the Japan Society for
Computational Engineering and Science, 14 (2009), pp. 2102-2123.
(7) He, X. and Luo, L.-S., "Lattice Boltzmann model for theincompressible Navier-Stokes equation," J. Stat. Phys., 88 (1997),
pp.927-944.
(8) Onishi, J., Chen, Y. and Ohashi, H., "Lattice Boltzmann simulationof natural convection in a square cavity," JSME Int. J. Ser. B-Fluids
Therm. Eng., 44 (2001), pp. 53-62.
(9) Saiki, Y. and Seta, T., "Numerical Analysis of Incompressible Flowsby the Lattice Boltzmann Method," Transactions of the Japan
Society of Mechanical Engineers Ser. B, 74 (2008), pp. 2124-2131.
(10) Inamuro, T., Ogata, T., Tajima, S., and Konishi, N., "A lattice
Boltzmann method for incompressible two-phase flows with large
density differences," J. Comput. Phys., 198 (2004), pp.628-644.
(11) Kajishima, T., Numerical Simulation of Turbulent Flows, Yokendo(1999), pp. 55-56.
(12) Takada, N. and Tsutahara, M., "Evolution of viscous flow around asuddenly rotating circular cylinder in the lattice Boltzmann
method," Comput. Fluids, 27 (1998), pp.807-828.

Fig. 11 The temperature distributions on horizontal cross section for
Pr=2,Qin= 0.15 l/min. and Tin= -10C at t = 0.848s.
５．結 論
本研究では，電子デバイス・機器が集積する各種装置や情報通信施設の冷却用消費電力の低減を目的として，デバイスに直接載せる小型冷却ユニット内の液単相流れと固体および流体内部の熱移動に関する数値シミュレーションを実施した．その結果，設定した流動条件（Re = 51.2～172.1，Pr = 1.0～4.0，Pe = 51.2～319.3）
で以下の結果を得た．
(1) 冷却速度（単位時間当たりの低下温度）は，ユニット温度（シミュレーションの固体温度）と冷却液温度（シミュレーショ
ンの流入流体温度）との差とともに増加する．
(2) 冷却速度は冷却液流量の増加によっても向上する．
(3) Prandtl 数 Pr が大きくなると冷却速度は悪くなるが，液流動開始後しばらくすると各Pr数条件で冷却速度は一定になる．
(4) ユニットは冷却液流路付近の上流側から先に冷却される．
(5) 流路周辺では流体-固体間の熱伝達により温度が低下し，より内側の固体部分から熱が熱伝導により流路に向かって移動して冷却が進展する．
(6) 流入境界で観察される最大圧力は流量とともに増加し，同じ流量では Re 数が低い（粘性が大きい）場合により大きくな
る．
上記結果から，今回想定した液体冷却ユニットは基礎的冷却性7

Copyright © 2011 by JSFM