數值模擬質子交換膜燃料電池的方法

韓麗輝,于春梅,馮根生

(北京科技大學 冶金與生態工程學院，北京 100083)

燃料電池是將燃料的化學能直接轉換成電能的裝置。燃料電池具有效率高、無污染、重量輕等優點，被稱為21世紀的綠色環保能源，具有很大的發展潛力和應用前景[1]。在燃料電池的研究中，實驗室的物理實驗和數值計算必不可少，可靠的計算程序不僅為物理實驗提供理論指導，還可以優化設計燃料電池、節約實驗成本、提高科研工作效率。數值模擬軟件Fluent在解決流體計算問題方面應用很廣泛，并且該軟件具有燃料電池模塊，但其質子交換膜燃料電池模塊是針對氫氧燃料電池開發的，當用戶的燃料電池是甲醇或其他質子交換膜燃料電池時，該模塊就不能完全適用，用戶必須使用UDF編寫自己的燃料電池程序。本文從Fluent數值求解機理、燃料電池問題基本設置、編寫用戶UDF、程序調試等方面，詳細闡述了如何利用Fluent軟件來數值求解質子交換膜燃料電池問題。

1 Fluent求解機理

1.1 基本數據結構

了解Fluent數據存取方式及程序求解機理是編寫UDF的前提和基礎。用戶自定義函數需要從Fluent求解器中讀取數據，數據的基本載體為網格。每個網格單元稱為cell，每個cell都有節點(node)、面(face)、網格中心點(cell center)，見圖1。每個網格都由面包圍而成。Thread 是一種數據結構，可以用來存儲一組網格和一組面的信息。與網格相關的數據類型為Node、face_t、cell_t 、Thread、Domain，通過這些數據類型定位到數據載體。本文中求解燃料電池問題的UDF程序會用到Fluent中已定義的一些基本數據宏[2-4]見表1，如網格中心點坐標，網格體積、流體密度、流體溫度、流體壓力、混合物組分、變量梯度、多孔介質孔隙度等。 表1中，ND ND為Fluent中全局變量，對于二維問題，該值為2，對于三維問題，該值為3；另外宏C_P(c,t)讀取到的壓力是表壓力(Gauge Pressure)，其絕對壓力等于該值加上操作壓力，Fluent中全局操作壓力變量為op_pres。

圖1 Fluent 中基本數據結構

表1 Fluent中部分基礎數據宏

1.2 方程求解過程

求解燃料電池問題需要求解自定義標量傳輸方程，Fluent中以壓力為基本項的耦合求解器求解過程如圖2所示。在迭代計算外要執行兩步運算：一是初始化，即用Fluent初始化界面中的數據來初始化求解變量，之后調用用戶自定義的初始化程序，后者數據會覆蓋前者數據；二是調用用戶自定義邊界條件程序(PROFILE UDF)。

求解迭代循環過程：首先從執行用戶定義的調整程序(ADJUST UDF)開始，之后耦合求解連續方程和動量方程；最后順序求解能量方程、組分傳輸方程、湍流方程、用戶自定義標量傳輸方程；所有方程求解結束后執行物料屬性更新程序(用戶自定義PROPERTY UDF)；最后檢查是否滿足收斂標準，不滿足進行下一次迭代計算，滿足則退出循環計算結束。

了解fluent控制方程求解過程對于用戶編寫自己的UDF非常重要，根據求解過程，用戶就能在合適的位置給自定義內存賦值。用戶感興趣的計算數據是通過自定義內存存儲的，自定義內存的賦值位置可以放在程序調整宏 、控制方程源項、物料屬性UDF，而不能放在迭代循環外邊的初始化及邊界條件的UDF 程序中。

圖2 Fluent中壓力基本項求解過程

2 燃料電池問題的基本設置

2.1 調入氫氧燃料電池模塊

通過文本框命令行調用質子交換膜燃料電池模塊。雖然Fluent提供的質子交換膜燃料電池為氫氧燃料電池，不完全適用甲醇燃料電池問題，但同屬于質子交換膜燃料電池，其共同之處是可以借用的，如燃料電池各層材料屬性、水合量擴散系數、水合量在催化層固定值等。

2.2 基本模型的選擇

自定義標量UDS及內存UDM。在界面中定義自定義標量UDS及自定義內存UDM個數，同時要保證界面中定義的數量大于UDF程序中需要的數量。在質子交換膜燃料電池中至少需要3個標量，即電子電勢、離子電勢及水含量。自定義內存個數根據用戶需要而定，至少包括體電流密度、電流密度矢量、水含量、電滲拖曳系數、反應熱等。

選擇層流模型、能量方程及組分傳輸模型。在組分傳輸模型中不考慮體化學反應，化學反應是通過標量及源項的設置體現出來。

2.3 流體屬性的設置

從Fluent物料數據庫中選擇甲醇空氣燃燒混合物,該混合物包括甲醇、氧氣、水、二氧化碳、氮氣。按順序設定流體組分，其索引號從0開始，依次為1、2、3、4等，該索引號在UDF中定義物料組分擴散系數時會被使用。流體屬性除了密度、比熱、黏度等基本參數外，還需要設定2個參數，且這2個參數對計算結果有很大的影響。一個參數是流體組分擴散系數。在多孔介質結構中，氣體分子往往被毛細管壁所阻礙，為了計算這些阻礙的影響，修正擴散通量，一般是通過修正擴散系數實現的，該系數主要影響組分傳輸方程。在質子交換膜燃料電池問題中計算組分擴散系數時可以使用稀釋近似模型[5]：

(1)

另一個參數是自定義標量擴散系數(UDS,diffusivity)。對于電子電勢和離子電勢來說，其擴散系數即在不同流體區域的電導率，交換膜的電導率與電池溫度及水含量關系[6]為

(2)

式中，σmem為交換膜的電導率(Ω-1·m-1)；β為膜的電導率系數；λ為膜中水含量；ω為膜的電導率指數。

擴散層、催化層的電導率為常數，集流板固相區域的電導率需要在固體材料屬性UDS diffusivity 中設定。

2.4 區域的設置

共有9個計算域，包括陰陽極集流板(2個)、陰陽極流體通道(2個)、陰陽極擴散層(2個)、陰陽極催化層(2個)、質子交換膜(1個)。其中集流板為固體域，流體通道、擴散層、催化層、膜為流體域。區域屬性的設置包括多孔介質、源項及固定值。膜、催化層和擴散層均為多孔介質區域，要選擇各自的多孔介質材料。多孔介質主要引起動量損失，區域為多孔介質域時，在動量方程中自動引入動量損失源項。計算該源項需要黏性損失系數、慣性損失系數及孔隙度。當流動是層流流動時，慣性損失可以忽略。燃料電池數值計算實際上是在計算流動問題基礎上加載電勢傳輸方程，同時由于電化學反應引起能量方程、組分方程、連續方程的源項發生變化，必須要加載源項(UDF 程序)。流體穿過多孔介質區域從陽極到陰極或從陰極到陽極，這種穿透運動有兩種驅動力，一是由于動量驅動的，二是由于濃度差驅動的，為了使兩極流體不相混合，消除動量驅動作用的影響，在膜域對流體運動的速度設為固定值0。另外膜的水含量影響膜的導電率及水的電滲效果，所以水含量只存在于催化層和膜，對于其他流體域水含量都需設為固定值0。

2.5 邊界條件的設置

陰陽極入口的質量流量、流體溫度，混合物中各組分質量比例、出口壓力、出口回流氣體溫度、集流板溫度、操作壓力都是可變量，可根據情況設定。陰陽極集流板電勢或電勢梯度如下：

陽極電子電勢?an：

?an=0

(3)

陰極電子電勢?cat：

?cat=Vcell

(4)

膜離子電勢?m：

(5)

式中n為膜厚。

2.6 求解控制及初始值的設置

壓力與速度耦合求解采用SIMPLE算法，采用固定循環模式，即按變量順序求解，使用BGSTAB方法增加求解過程的穩定性。給自定義標量電子電勢和離子電勢賦予合理的初始值，是計算是否順利進行的關鍵，如果電勢初始值不合理，自定義標量殘差曲線出現振蕩或計算出現錯誤而停止。合理的電勢初始值能夠使得陰陽極都有合理的過電勢，即電子和離子運動的驅動力，這樣才能使電化學反應連續發生。一般離子電勢的初始值為負數，如-0.1。

3 編寫燃料電池程序

在Fluent軟件提供的動量、能量、連續、組分傳輸等方程的基礎上加載電子電勢及離子電勢標量傳輸方程，并對原有方程加載源項，并修改混合流體的相關屬性，最后數值求解方程組，得到流場、質量組分、電流密度、溫度場的分布情況。

3.1 質子交換膜燃料電池基礎

質子交換膜燃料電池(PEMFC)結構如圖3所示，其陰陽兩極各包括集流板(collector)、流道(gas channal)、擴散層(diffusion layer)和催化層(catalyst layer)，陰陽兩極中間是電解質膜(electrolyte membrane)[7-8]。燃料從陽極流道入口流入，氧化劑從陰極流道入口流入，陰陽兩極入口分別在燃料電池的兩個端部，流體相對而行。擴散層和催化層為多孔介質區域，流體擴散到催化層后發生電化學反應，陽極電化學反應產生的氫離子通過電解質膜進入陰極，電子則通過陽極擴散層及集流板和外電路進入到陰極集流板，再經陰極擴散層，進入陰極催化層，在那里和氧化物及氫離子發生陰極電化學反應。陰陽兩極不斷供應燃料和氧化劑，電化學反應就不斷的進行下去，外電路中的負載就得到了電能的供應，這就是簡單的燃料電池將燃料化學能轉化為電能的過程。

甲醇燃料電池的陽極、陰極催化層電化學反應方程式分別為：

陽極催化層：

CH3OH+H2O=CO2+6e-+6H+

(6)

圖3 PEMFC 結構圖

陰極催化層：

(7)

陽極催化層電化學反應產生二氧化碳、氫離子和電子，二氧化碳從陽極出口流出，電子從外電路到達陰極，氫離子通過電解質膜到達陰極，這樣就形成一個閉合回路[9]。陰極氧氣和電子及氫離子反應產生水，水會反向擴散到陽極。在質子交換膜燃料電池中，氫離子的傳輸離不開水，只有膜的納米尺度空隙中存在足夠的水時，水合氫離子才能夠從陽極到達陰極。攜帶氫離子的水為離子拖曳水，這種現象為電滲拖曳[10]。在甲醇燃料電池中，不僅有水的電滲拖曳，還有甲醇的電滲拖曳和甲醇的擴散，到達陰極的甲醇和陰極的氧反應，產生寄生電流。

3.2 電化學模型

在PEMFC的陰陽極催化層不斷發生電化學反應，不斷產生和消耗電子和氫離子，反應發生的驅動力是同一網格上固相電勢(電子電勢)與膜電勢(離子電勢)之間的電勢差，即活化過電勢。電子電勢傳輸主要在固相區域中，包括集流板、擴散層及催化層，其傳輸方程見式(8)。離子電勢的傳輸發生在催化層及膜區域中，其傳輸方程見式(9)。電勢傳輸方程中的擴散系數即電導率，包括電子在催化層、擴散層、集流板中的電導率，及離子在催化層和膜中的電導率。傳輸方程中的源項R為體電流密度，即單位時間單位體積產生的電量。體電流的大小與反應物濃度、活化過電勢、系統溫度、離子交換膜厚度等有關[11-12]。

·(σsol?sol)+Rsol=0

(8)

(9)

式中，σ為電導率(Ω-1·m-1)，下標sol表示固體(陰陽電極)，m代表膜；?為電勢(V)；R為體電流密度(A·m-3)。

采用簡化的Butler-Volmer 公式(即Tafel 公式)來計算陽極和陰極的交換電流密度分別為：

(10)

(11)

式中，jref為參考交換電流密度(A·m-2)；ζ為反應區比表面積(1/m)；[A]和[A]ref分別為反應物濃度及參考濃度 (mol·m-3)；γ為濃度指數；ηan為過電勢；α為傳遞系數；F為法拉第常數(96 485 C/mol)。

在陽極催化層，電子體電流密度等于負的陽極交換電流密度，離子體電流密度等于陽極交換電流密度；在陰極催化層，電子體電流密度等于陰極交換電流密度，離子體電流密度等于負的陰極交換電流密度。陽極陰極的活化過電勢不同，有：

陽極催化層：

Rsol=-Ran(<0)，Rm=+Ran(>0)

陰極催化層：

Rsol=+Rcat(>0) ，Rm=-Rcat(<0)

陽極過電勢：

ηan=?sol-?m

(12)

陰極過電勢：

ηcat=?sol-?m-Vopen

(13)

3.3 質量源項

流體混合物組分發生變化的地方，需要在連續方程和組分傳輸方程中增加源項。燃料電池陽極催化層和陰極催化層分別發生電化學反應，從反應方程式(6)和(7)可以看出，陽極催化層消耗掉甲醇和水，產生二氧化碳，陰極催化層消耗掉氧氣，產生水，當考慮寄生電流時，陰極還會有二氧化碳產生。陽極催化層消耗的甲醇包括電化學反應消耗、電滲拖曳和甲醇擴散，其質量源項表達式見式(14)；陽極水的變化包括電化學反應和電滲拖曳消耗的水及反向擴散來的水，其質量源項表達式見式(15)；陽極電化學反應產生二氧化碳，其質量源項表達式見式(16)。陰極催化層消耗的氧氣包括產生主電流需要的氧氣和產生寄生電流需要的氧氣，其質量源項表達式見式(17)；陰極水的變化包括主電流反應產生的水、寄生電流反應產生的水、電滲拖曳作用帶來的水和反向擴散減少的水，其質量源項表達式見式(18)；陰極產生的二氧化碳由于寄生電流而存在，其表達式見式(19)。在連續方程(mass 方程)中加入變化組分的質量源項。

Dch3ohρ(2Y))

(14)

(15)

(16)

(17)

(18)

(19)

式中，S為質量源項(kg·m-3·s-1)；M為分子量(kg·m-3) ；ρ為混合物密度(kg·m-3) ；Y為組分質量分數；Rp為寄生電流密度(A·m-3)；ndr為電滲拖曳系數；Dλ為水含量擴散系數(mol·m-1·s-1) ；Dch3oh為甲醇擴散系數(m2·s-1) 。

3.4 能量源項

由于電流的產生，使得燃料電池除流道外的各個區域的能量方程中都需要增加源項。在催化層，電化學反應的焓變并沒有完全轉化成電能，為了克服活化能壘有一部分能量損失掉了，催化層能量方程源項見式(20)，在集流板、擴散層和交換膜中，能量方程源項為導電介質的散熱量，見式(21)，式中電流密度為單位時間通過單位面積上的電荷，該量為矢量，在某一方向上分量大小為該方向上電勢梯度與電導率的乘積。如果考慮到水蒸氣發生相變，那么要在擴散層和催化層的能量源項中增加水的汽化潛熱。

Sh=-Δhrect-Ran,catηan,cat

(20)

Sh=I2Rohm

(21)

式中，Δhrect為電化學反應焓變(W·m-3)；I為電流密度(A·m-2)，I=σ?；Rohm為電阻，

3.5 自定義內存

在不同自定義函數之間共用的變量值存儲在自定義內存中，因此給自定義內存賦值位置在求解燃料電池問題中極為中重要，要根據求解過程中不同宏的調用順序，在該內存值在被引用之前為其賦值。在燃料電池問題中，至少需要如下自定義內存：過電勢、交換體電流密度、水含量、膜的電導率、水的電滲拖曳系數、電流密度及其在x、y、z方向上的分量，以及在計算過程及結果中感興趣的變量等。過電勢在計算交換體電流密度時使用，在電勢傳輸方程源項中賦值；交換體電流密度在各個方程的源項中都可能用到，所以定義在adjust宏中賦值；水含量在計算膜電導率、水含量擴散系數及水的電滲拖曳系數時都會用到，所以在混合物屬性宏中賦值；膜的電導率及水的電滲拖曳系數也在屬性宏中賦值；電流密度及其分量只在能量方程源項中引用，且除流道域外各個區域的能量方程均勻源項，所以在能量方程源項中給電流密度賦值即可。在計算過程陽極陰極的電流不斷接近，收斂時達到相等，我們既關心電流過程值便于調試程序，也關心計算結果電流的大小，所以在自定義內存中要有電流，電流是催化層所有網格體電流密度的總和，在adjust宏中給電流賦值。

4 程序調試

電化學反應主要通過電勢標量傳輸方程及其他方程增加源項體現出來的。程序調試時，首先不考慮電化學反應，即不求解電勢標量傳輸方程也不增加源項的情況下，確保混合物流場的正確性。要保證陽極入口流入的流體在陽極流道、陰極入口流入的流體在陰極流道，一是要設定膜域中由于壓力差引起的宏觀流動的速度為零，二是設定合理的混合物擴散系數。待無電化學反應的流場計算正確后，再加載adjust 宏及電勢傳輸方程中源項及質量源項，調試電流密度場，最后加載能量方程中源項，調試溫度場。

5 計算結果

開路電壓為0.95 V，進料溫度353 K，陰陽極參考電流密度為1 000 A/m2，操作壓力為101 325 Pa，陽極入口流體質量流量為5×10-6kg/s，陽極入口流體質量流量為6×10-7kg/s，陰陽極電荷傳遞系數為2，電池擴散層、催化層及膜的氣孔率為0.5，陽極濃度指數為0.5，陰極濃度指數為1。流道為單通道時計算得到的電流密度場見圖4、溫度場見圖5、甲醇質量場見圖6、水分質量場見圖7。

圖4 電流密度分布

圖5 溫度分布

圖6 甲醇質量組分分布

圖7 水分質量組分分布

6 總結

用Fluent 流體計算軟件求解質子交換膜燃料電池問題，關鍵在于設定正確的自定義電勢標量擴散系數(即電導率)、標量傳輸方程源項(體電流密度)、連續方程及組分傳輸方程中質量源項、能量方程中能量源項，還要設定合理的自定義內存變量及合適的賦值位置。

[1] 奧黑爾.燃料電池基礎[M].王曉紅，黃宏，譯.北京：電子工業出版社，2007.

[2] 韓占鐘，王敬，蘭小平.Fluent流體工程仿真計算實例與應用[M].北京：北京理工大學出版社，2004.

[3] 王瑞金，張凱，王剛.Fluent技術基礎與應用實例[M].北京：清華大學出版社，2007.

[4] 張建文，楊振亞，張政.流體流動與傳熱過程的數值模擬基礎與應用[M].北京：化學工業出版社，2009.

[5] Sukkee U m,Wang C Y,Chen K S.Computational Fluid Dynamics Modelingof Proton Exchange Membrane Fuel Cells[J].J Electrochemical Society,2000，147(12):4485-4493.

[6] Springer T E,Zawodzinski T A,Gottesfeld S.Polymer Electolyte Fuel CellModel[J].J Electrochemical Soc,1991，138(8):2334-2342.

[7] 熊濟時，肖金生，潘牧.質子交換膜燃料電池流場模擬與結構尺寸優化[J].武漢理工大學學報：交通科學與工程版,2009,33(6)：534-537.

[8] 馬利軍，林才順，薛方勤，等.質子交換膜燃料電池流場截面設計及三維模擬[J].濕法冶金，2008,27(1):52-55.

[9] 劉桂成，張浩，王一拓，等.液體進料直接甲醇燃料電池工況的匹配及優化[J].電池，2012,42(1):7-10.

[10] 潘科，齊亮，劉洋，等.質子交換膜的電滲系數[J].電池，2008,38(2):75-78.

[11] Kulikovsky A A,Divisek J,Kornyshev A A.Modeling the Cathode Com-partment of Polymer Electrolyte Fuel Cells: Dead and Active Reaction Zones[J].J Electrochemical Society,1999，146(11):3981-3991.

[12] Mazumder S，Cole J V.Rigorous 3-D Mathematical Modeling of PEM Fuel Cells II.Model Predictions with Liquid Water Transport[J].J Electrochemical Soc,2003，150(11):A1510-A1517.