應用晶格波茲曼法與場協同理論於增強對流熱傳之研究

行政院國家科學委員會補助專題研究計畫期末報告

計畫名稱：應用晶格波茲曼法與場協同理論於增強對流熱傳之研究

計畫類別：個別型計畫

計畫編號：NSC96-2221-E006-168-MY3

執行期間：九十六年八月一日至九十九年七月三十一日 計畫主持人：陳朝光教授 國立成功大學機械工程學系 執行單位：國立成功大學機械工程學系

中華民國九十九年十月二十五日

行政院國家科學委員會補助專題研究計畫期末報告

計畫名稱：應用晶格波茲曼法與場協同理論於增強對流熱傳之研究 計畫編號：NSC96-2221-E006-168-MY3

執行期間：九十六年八月一日至九十九年七月三十一日 計畫主持人：陳朝光教授 執行單位：國立成功大學機械工程學系

摘要

晶格波茲曼法(Lattice Boltzmann Method, LBM)近年來被應用在 各種物理領域的模擬與分析越來越多，尤其是在計算流體力學方面。

本計劃即利用晶格波茲曼法模擬質子交換膜燃料電池雙極板之流 道，模擬範圍包含了流體在雙極板流道和擴散層中流動、擴散與熱傳 現象，屬於低雷諾數二維穩態不可壓縮的流場。並藉由實驗結果來驗 證模擬之正確性。將置入障礙物後的流道視為波浪型流道且流場之熱 傳增強分析利用場協同理論來驗證得到的數值解趨勢，透過速度場與 溫度梯度場之間的協同角越小，表示其協同程度越好，進而幫助熱傳 量的增加。

關鍵字：晶格波茲曼法、場協同理論、質子交換膜燃料電池、障礙物

ABSTRACT

This study applies the Lattice Boltzmann Method (LBM) to simulate gas flow characteristic in bipolar plate channel of Proton Exchange Membrane Fuel Cells (PEMFCs). The calculated domain includes fluid flow, diffusion, and heat transfer performance in bipolar plate and the porous layer which is been considered as a two-dimensional incompressible steady flow under low Reynolds number. And compare the experimental results with the simulation results. Furthermore, the present numerical results are consistent with the field synergy principle, which states that the convective heat transfer is enhanced when the intersection angle between the velocity and the temperature gradient vector decreased by inserting obstacles to fluid field. The results of heat transfer improved significantly.

Keyword : Lattice Boltzmann Method, field synergy principle, PEMFC, Obstacles

前言

當全球社會正面臨空前的能源危機與環境污染問題時，許多具前 瞻性的國家早已研發高效率且低污染的新能源，在眾多替代能源技術 裡，燃料電池技術脫穎而出，成為全球最熱門的研究方向。燃料電池 是結合燃料與氧化劑而產生電能的一種電化學裝置，其典型的燃料為 氫氣而氧化劑為氧氣或空氣，操作時的反應產物為水和熱。在氫氣轉 換成電能時沒有燃燒作用，整個發電過程是高效率、無污染且安靜的 系統。當前質子交換膜燃料電池(PEMFC)相關研究大致可區分為流道 設計、導電材料、質子交換膜性能改良與催化反應機制等研究方向，

其中流道設計多著重於如何增加在雙極板內燃料可流動之區域，以期 在較長的流動時間內，增加流入擴散層中之燃料。但在流道內的流動 過程中，燃料進入擴散層內主要的機制為擴散作用，其流動速度極為 緩慢，若能在流動之中增加對流機制的效應，可有效提升燃料傳輸至 擴散層的量與速度。本研究就專以晶格波茲曼法模擬質子交換膜燃料 電池雙極板流道內的熱質傳現象，模擬範圍包含了雙極板流道與擴散 層，並設計不同的雙極板流道形狀，期望利用流道形狀的改變，除了 本身的擴散效應外，加入流體流動時的對流效應，增強流體流入擴散 層內的量與速度，例如放置半圓形障礙物在壁面(可視為波浪型流道) 取代直線型壁面。針對不同類型壁面對於流體的流動情形之影響與熱

傳現象做詳細的討論，經由此數值分析模擬的方法不僅可加速開發新 ㄧ代的流道，也可對燃料電池內的微觀物理現象能更進ㄧ步的理解，

這對於從事燃料電池研究工作而言，是一個極有力且不可或缺的工 具。

晶格波茲曼法(Lattice Boltzmann Method)是近年來發展的一種計 算流體力學方法，與一般以流體運動的統御方程式為出發點的計算流 體力學方法不同，一般計算流體力學方法是直接求解Navier-Stokes 方程式中的巨觀物理量，而晶格波茲曼法則是以流體的微觀模型為出 發點，是一種全新的計算流體力學的方法。

基本模型

在LBM中，最常使用的是Qian於1992年所提出的DnQb系列模 型，其中n代表空間維數，b代表離散速度數目。本研究採用最普遍使 用的D2Q9模型來模擬二維流場，D2Q9模型由正方形晶格組成（如圖 一所示），在每個正方形格點上都有三種型態的質點：靜止的質點

（

i

 ）、沿著水平或垂直方向移動的質點（0

i

1, 2,3, 4）、沿著對角 線方向（與水平或垂直方向成45度角）的質點（

i

5,6,7,8）。

晶格速度向量



e

i的定義如下

 

^0,0

 

e

i ，

i

 0

 ¹   ¹ 

cos ,sin

2 2

 

 

 

  

 



i

i i

e

，

i

1, 2,3, 4

 ¹   ¹ 

2 cos ,sin

2 4 2 4

   

       

        

   

 



i

i i

e

，

i

5,6,7,8

圖一

D2Q9

模型晶格速度向量示意圖

晶格波茲曼法的演化方程式

:



^{( , ) ( , )}



) 1 , ( ) ,

(x e t t t f x t f x t f⁽⁰⁾ x t

f_i _i       _i 









^{( , ) ( , )}



) 1 , ( ) ,

(x e t t t g x t g x t g⁽⁰⁾ x t

g_i _i       _i 







其中 f_i是粒子分佈函數，g_i為內能密度分佈函數，x是粒子的位置向 量，t為時間，_代表密度分佈函數趨近於局部平衡態 f⁽⁰⁾的弛豫時間

(relaxation time)

。

平衡態密度分佈函數為 ⁽⁰⁾  _  ²_^

2 1 3 9 ) 4 ,

(x t u

f_i    ，

i  0



     

 ^{2 2}

) 0 (

2 ) 3 2( ) 9 ( 3 9 1

) 1 ,

(x t e u e u u

f_i        ，

i  1, 2,3, 4



     

 ^{2 2}

) 0 (

2 ) 3 2( ) 9 ( 3 36 1

) 1 ,

(x t e u e u u

f_i        ，

i  5,6,7,8

平衡態內能密度分佈函數為

 0 ²

2

2

i 3

g u

c

 





，

i  0

 0

   

^{2 2}

2 4 2

3 3 9 3

9 2 2 2 2

i i

i

e u e u u

g c c c

 ^   ^

     

 

 

    

，

i

1, 2,3, 4

 0

   

^{2 2}

2 4 2

9 3

36 3 6 2 2

i i

i

e u e u u

g c c c

 ^   ^

     

 

 

    

，

i

5,6,7,8

巨觀物理量可表示如下





i

fi



^



i i if e

V 



i i

   g

多孔性模型

考慮多孔性質流的統御方程式時，在動量方程式裡會多出線性的

Darcy

項與非線性的

Forcheimer’s

項，因此我們可以將統御方程式裡

的外力項F_i加入晶格波茲曼方程式



_,

  

_, i

 

^, i^{ eq}

 

^,

i i i i

f x t f x t

f x e t t   t f x t  tF



        

  



²



2 4

1 F F : 1 2

i i s

i

i i

s s

u e e c I F e

c c

   

   

 

     

  

  



4 / 9 0 1/ 9 1,2,3,4 1/ 36 5,6,7,8

i

i i i



 

   

 



當加入外力項時，動量守恆式則變為 F

i i 2

i

u e f  t



^



^ 



^，因為

^F

^包含

了速度的非線性平方項，因此為了求解速度

u

，需先求解假定速度

v

，

而

v

可定義為 G

i i 2

i

v e f  t



^



^ 



^{，最後在套上}

^u

^與

^v

^{之間的關係式。}

2

0 0 1

u v

c c c v

  

 

 ，其中 ₀ 1

2 1 2

c t

K

  

 

   、 ₁

2 F c t

K









。

數值模擬結果與討論 流場速度邊界條件：

(1)

入口邊界

在本研究的流場中，入口之速度分佈為一均勻

(uniform)

分佈 U U in ，V 0

(2)

壁面邊界

壁面皆為無滑移邊界 U0，V 0

(3)

出口邊界

出口為固定壓力邊界且速度採取外插方式得到 溫度邊界設定：

(1)

入口邊界

入口邊界採取等溫邊界T T _in

(2)

下壁面邊界

下壁面邊界為等溫邊界T T _w

(3)

出口邊界

出口邊界採取熱全展邊界T_i T_i1 速度場與溫度場的收斂標準：

   

 

, ,

, 6

, ,

, ,

10 ,



_

 







   

 

i j i j

i j

i j i j

u x t u x t

u x t

   

 

, ,

, 6

, ,

, ,

10 ,



_

 







 



i j i j

i j

i j i j

T x t T x t

T x t

本節所要討論的是在流道下壁面放置半圓形障礙物

(

可視為波浪狀流 道

)

，加入障礙物的主要目的是想藉由障礙物阻隔而產生強制對流作 用，一開始流體在流道上游是以擴散作用傳輸至擴散層，當流體流經 障礙物位置時，遭受障礙物阻隔無法繼續順沿流道方向流動，迫使流 體經由強制對流作用進入擴散層參與反應。

首先圖二為具半圓形障礙物之雙極板流道示意圖，圖中每一個半 圓形障礙物的半徑皆為

H/6

，兩兩圓心之間的距離為

3H

，第一個半 圓形障礙物的圓心距離入口處為

20H/3

，最後一個半圓形障礙物圓心 距離出口處為

H

。圖三為

Re=80

的條件下具半圓形障礙物之流道內速 度向量分佈的局部放大圖，當流體通過半圓形障礙物時，流道忽然變 窄，加速了流體流動速度，在半圓形障礙物的後方出現一個小迴流

區，此迴流區的大小會伴隨著

Re

數減少而跟著縮小。圖四為

Re=80

的條件下具半圓形障礙物之流道內流線的局部放大圖，可以看出流線 經過半圓形障礙物的上方時受到擠壓，加速流體流動。圖五為

Re=80

的條件下具半圓形障礙物之流道內無因次溫度分佈圖，障礙物部份所 設定的溫度與下壁面溫度一樣，藉由下壁面凸起的障礙物可以幫助流 場的熱量由此傳出。圖六為具半圓形障礙物之雙極板流道內

Re

數對 壁面

Nu

數分佈之影響，由圖六可看到壁面

Nu

數有五個週期性變化 的趨勢，這個現象與流場的下壁面幾何形狀接觸位置有關，當流體流 經過半圓形障礙物時截面忽然變窄，因此加速了流體速度產生

y

方向 的速度，進而壓縮了熱邊界層，使得局部

Nu

數提高，並在障礙物的 最高點達到最大值，流體通過最窄截面後流道漸漸擴張，而

Nu

數也 跟著下降。在

Re=80

的條件下多加了半圓形障礙物之後的壁面

Nu

數 的最大峰值相較於無障礙物的流場大約增加了

196%

。

以下利用場協同理論來證實數值模擬的結果趨勢，根據場協同理 論若要增加對流熱傳效果則無因次積分值應該盡可能大，意即速度向 量與溫度梯度向量夾角應盡可能的小。圖七為流道內有無半圓形障礙 物情況下雷諾數對平均協同角的影響，協同角越小代表流場中速度場 與溫度梯度場之間協同程度越好，障礙物在流場中起了擾動的效果，

因而改善速度場與溫度梯度場之前的協同程度，具半圓形障礙物之流

道其平均協同角在不同雷諾數均小於無障礙物的流道，在

Re=80

的 條件下兩者之間的平均協同角大約相差

2

度，由此可證實模擬結果的 正確性。

圖八為實體波浪型雙極板之上視圖，其成分是不銹鋼合金纖維。

圖九為質子交換膜燃料電池雙極板之整體組合，測試機台採用

Bean

公司型號

300M

設備，如圖十所示。測試機台主要包含氣體供應系統、

流量控制系統、溫度控制系統、增濕系統與電子負載系統等五部份。

實驗之氣體供應系統主要供應氧氣、氫氣與空氣作為燃料電池陰極與 陽極反應氣體，而氮氣對燃料電池而言，其主要功能在於實驗前、後 清除管路中所殘留氣體與雜質。燃料氣體與氧化劑採用高純度氫氣與 氧氣，其純度為

99.99 %

以上。實驗中流量控制系統主要包含以化學 計量數

(Stoichiometry)

與最小流量

(Minimum Flow)

方式控制氣體進口 流量；在溫度控制系統方面，係使用加熱棒、

T

型熱電偶

(T-type

Thermocouple)

與溫控器等三部分來控制燃料電池內部工作溫度；在增

濕系統方面，係使用外部潤濕法中的水瓶氣泡對氣體燃料做加濕工 作，並配合溫度控制器來調節加濕瓶內部水溫，以確保氣體通過加濕 瓶達到所需入口氣體加濕溫度，以及在此溫度下氣體所含飽和蒸汽 壓；在電子負載系統方面，最大輸出功率為

300 W

，最大電流為

120

A

，最大電壓為

60 V

。

圖二 具半圓形障礙物之雙極板流道示意圖

x/H

5.5 6 6.5 7 7.5

圖三 具半圓形障礙物之雙極板流道內局部放大速度向量分佈圖

(Re=80)

x/H

5.5 6 6.5 7 7.5

圖四 具半圓形障礙物之雙極板流道內局部放大流線圖

(Re=80)

x/H

0 5 10 15 20

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1

圖五 具半圓形障礙物之雙極板流道內無因次溫度分佈圖

(Re=80)

0 4 8 12 16 20 x/H

0 2 4 6

Nu

Re=80 Re=70 Re=60 Re=50 Re=40

圖六 具半圓形障礙物之雙極板流道內

Re

數對下壁面

Nu

數分佈之 影響

40 50 60 70 80 Re

83 84 85 86 87

intersection angle, degree

type1 type2

圖七

Re

數對平均協同角之影響

(type1:

平板流道

; type2:

波浪型流道

)

圖八 不銹鋼合金纖維製作完成之雙極板實體圖

圖九 單一燃料電池實體圖

圖十 質子交換膜燃料電池量測設備圖

參考文獻

[1] Liu, H.C., Yan, W.M., Soong, C.Y., Chen, F., “Effects of baffleblocked flow channel on reactant transport and cell performance of a proton exchange membrane fuel cell,” Journal of Power Sources 2005, 142, 125-133.

[2] Soong, C.Y., Yan, W.M., Tseng, C.Y., Liu, H.C., Chen, F., Chu, H.S., “Analysis of reactant gas transport in a PEM fuel cell with partially blocked fuel flow channels,”

Journal of Power Sources 2005, 143, 36-47.

[3] Yan, W.M., Chen, F., Wu, H.Y., Soong, C.Y., Chu, H.S., “Analysis of thermal and water management with temperature-dependent diffusion effects in membrane of proton exchange membrane fuel cells,” Journal of Power Sources 2004, 129, 127-137.

[4] Guo, Zhaoli, & Zhao, T. S., “Lattice Boltzmann model for incompressible flows through porous media”, Phys. Rev. E, Vol. 66, pp. 036304, 2002.

[5] Guo, Zhaoli, & Zheng, Chuguang, “An extrapolation method for boundary conditions in lattice Boltzmann method”, Phys. Fluids, Vol. 14(6), pp. 2007-2010, 2002.

[6] Guo, Z. Y., & Tao, W. Q., & Shah, R. K., “The field synergy (coordination) principle and its application in enhancing single phase convective heat transfer”, Int. J.

Heat Mass Transfer, Vol. 48, pp. 1797-1807, 2005.

[7] Gupta, Abhishek K., & Sharma, Atul, & Chhabra, Rajendra P., & Eswaran, Vinayak, “Two-dimensional steady flow of a power-law fluid past a square cylinder in a plane channel: Momentum and heat-transfer characteristics”, Ind. Eng. Chem. Res., Vol. 42, pp. 5674-5686, 2003.

[8] He, Xiaoyi, & Chen, Shiyi, & Doolen, Gary D., “A novel thermal model for the lattice Boltzmann in incompressible limit”, J. Comput. Phys., Vol. 146, pp. 282-300, 1998.

[9] Lim, C.Y., & Shu, C., & Niu, X. D., & Chew, Y.T., “Application of lattice Boltzmann model to simulate microchannel flows”, Physics of Fluids, Vol. 14(7), pp.

2299-2308, 2002.

[10] McNamara, Guy R., & Zanetti, Gianluigi, “Use of the Boltzmann equation to simulate lattice-gas automata”, Physical Review Letters, Vol. 61(20), pp. 2332-2335, 1988.

[11] Mei, Renwei, & Luo, Li-Shi, “An Accurate curved boundary treatment in the lattice Boltzmann method”, J. Comput. Phys., Vol. 155, pp. 307-330, 1999.

[12] Munson, Bruce R., & Young, Donald F., & Okiishi, Theodore H., “Fundamentals of fluid mechanics”, Wiley, New York, 1998.

[13] Noble, David R., & Chen, Shiyi, & Georgiadis, John G.., “A consistent hydrodynamic boundary condition for the lattice Boltzmann method”, Phys. Fluids, Vol. 7(1), pp. 203-209, 1995.