廖梓豪，魏 琳，索振邦，石偉玉，蔣方明

質子交換膜燃料電池冷啟動時陰極催化層內多場耦合的介孔尺度模型*

廖梓豪1,2,3,4，魏 琳1,2,3，索振邦5，石偉玉5，蔣方明1,2,3?

（1. 中國科學院廣州能源研究所，廣州 510640；2. 中國科學院可再生能源重點實驗室，廣州 510640； 3. 廣東省新能源和可再生能源研究開發(fā)與應用重點實驗室，廣州 510640；4. 中國科學院大學，北京 100049； 5. 上海捷氫科技有限公司，上海 201800）

為詳細解析質子交換膜燃料電池（PEMFC）零下溫度啟動過程，建立了電池冷啟動時多場耦合過程的介孔尺度數(shù)值模型。幾何模型基于隨機網(wǎng)格法（SGM）重建的催化層介孔結構，數(shù)學模型描述了物質傳輸、電荷傳輸、電化學反應和汽?冰相變過程。數(shù)值分析了電池冷啟動過程中催化層內冰的生成和演化，重點探討了冰的生成及形貌對電池性能的影響，推導出電化學活性反應面積與結冰量的關系式。

燃料電池；數(shù)值模擬；介尺度；催化層；冷啟動

0 引 言

質子交換膜燃料電池（proton exchange membranefuel cell, PEMFC）具有能量轉化效率高、啟動快、零排放等優(yōu)點，應用前景廣闊[1]。然而，PEMFC的規(guī)?；褂萌匀幻媾R諸多挑戰(zhàn)。例如，電池自零下溫度環(huán)境啟動時，可能會因結冰導致反應位點被覆蓋或者反應（產(chǎn)）物傳輸通道被堵塞而無法啟動。

數(shù)值模型是研究PEMFC內多場耦合過程的重要手段和工具。MAO等[2]首先建立了PEMFC的三維冷啟動模型，研究了冷啟動時電流密度、含水量和冰的三維分布特征；JIANG等[3]為了解PEMFC冷啟動時非等溫特性，研究了因溫度升高導致電?熱?氣（水）傳輸物性的改變及其對電池冷啟動性能的影響。盡管宏觀尺度冷啟動模型[3-5]已經(jīng)實現(xiàn)PEMFC冷啟動過程的仿真，但是尚不能夠對結冰形貌等進行動態(tài)可視化，宏觀模型中催化層（catalyst layer, CL）內含冰量對物質輸運以及電化學有效反應面積的影響關系式都是基于純粹的經(jīng)驗關系式，缺乏必要的微觀基礎。相對于宏觀均質模型，介孔模型基于電極的介孔微結構建模[6-7]，模型結果可揭示電極的構效關系，可用于設計或優(yōu)化電極微孔結構，還可獲得更準確的有效傳輸物性等宏觀參數(shù)（或關系式）[8]。WANG等[9]首次利用理想規(guī)則結構的催化層建立了PEMFC過程介孔尺度數(shù)值模型，并研究了催化層形貌（孔徑、厚度和孔隙率）對電池性能的影響，隨后模型進一步發(fā)展為基于隨機催化層微結構[10]。

微孔結構乃孔尺度模型的基礎。微結構重建方法主要分為實驗重建法[11]和數(shù)值重建法[6,12]。與實驗重建法相比，數(shù)值重建法具有節(jié)省時間和降低成本等優(yōu)點。隨機網(wǎng)格法[13]是數(shù)值重建法的一種，作為Monte Carlo模擬在Ising模型的一種擴展，能利用電極的統(tǒng)計學數(shù)據(jù)對電極結構進行重建，計算速度較快，重建結構能較好地重現(xiàn)原結構的統(tǒng)計特征。

本文基于隨機網(wǎng)格法重建催化層介孔結構，針對PEMFC冷啟動時催化層內發(fā)生的電?熱?氣（水）傳輸、電化學反應和結冰的多場耦合過程，建立了介孔尺度模型。研究了冷啟動過程中，冰的形成和分布以及冰的形成對傳質及電化學反應的影響，詳細分析了PEMFC冷啟動行為特性。

1 陰極催化層微結構重建

擬重建的催化層厚度為15 μm，由固體混合相（mixed solid, MS）（包括炭黑、Pt和離子聚合物）和孔相（pore, P）組成，各相體積分數(shù)見表1。采用隨機網(wǎng)格法[14]（stochastic grid method, SGM）重構微結構，重建結構的兩點自相關函數(shù)由實際催化層的SEM圖（圖1）提取，SEM圖的像素位點尺寸為31.716 nm，因此在厚度方向（）上設置473個數(shù)值單元來代表厚度為15 μm的催化層；在另外兩個方向（、）上各設置20個數(shù)值單元。重建結構如圖2，其中紅色代表固體混合相、淡藍色為孔相。

表1 催化層的各相組成

由隨機行走模型[15]計算得到重建結構的扭曲率為1.505，與實際測試結果大致吻合[16]。采用Hoshen-Kopelman聚類分析算法對重構微結構進行連通性分析，連通的孔相和MS相分別占各自總量的97.26%和99.39%。

圖1 陰極催化層的SEM圖

圖 2 重建結果（紅色為固體混合相，淡藍色為孔相）

圖3所示為重構微結構的各相體積分數(shù)（截面平均）沿電極厚度方向的變化曲線，各相體積分數(shù)在平均體積份額附近波動，波動幅度小于10%。圖4將重構微結構的孔徑分布與壓汞實驗（mercury injection experiment, MIP）結果進行了對比。實驗曲線中的兩個峰分別代表燃料電池催化層內的微孔和大孔，重建結構中的孔徑分布接近目標電極的大孔區(qū)域；本文重建結構的解析度小于實際的微孔尺寸，故重建結構無法呈現(xiàn)微孔。

圖3 各相體積分數(shù)沿電極厚度方向的變化

圖4 重建催化層的孔徑分布與實驗結果的對比

2 數(shù)值模型

2.1 模型假設

陰極催化層位于質子交換膜（proton exchange membrane, PEM）和氣體擴散層（gas diffusion layer, GDL）之間（如圖5示），= 0處為催化層與氣體擴散層交界面，=L（L= 15 μm）處為催化層與PEM交界面；GDL為氧氣、水蒸氣和電子傳輸提供通道，PEM為質子和膜態(tài)水傳輸提供輸運路徑。從GDL側進入的氧氣與陽極反應產(chǎn)生的電子（經(jīng)外電路后由GDL側進入）和質子（通過PEM側進入）在催化層內相遇，在催化劑作用下在活性位點（相界面）發(fā)生電化學反應。產(chǎn)物水一部分以水蒸氣形式向GDL側擴散；另一部分以膜態(tài)水的形式向PEM側擴散（回擴散）。此外，電滲作用會使陽極及PEM內部分膜態(tài)水遷移到陰極催化層內。當催化層內的水（水蒸氣）出現(xiàn)局部過飽和時，將發(fā)生相變而形成冰。

模型主要假設如下：①結冰為汽?固凝華過程，且結冰只發(fā)生在孔隙之中；②電化學反應發(fā)生在孔相和固體混合相之間的交界面；③為專注于研究水的輸運對冷啟動的影響，模型僅考慮等溫冷啟動過程[16-17]；④忽略催化層內超低速流動對流傳質；⑤忽略催化層中微孔對傳質與電化學反應的貢獻。

圖5 陰極催化層模型的電荷及物質輸運示意圖

2.2 控制方程

模型求解5個控制方程，分別是氧氣濃度、水蒸氣濃度、膜含水量方程及質子、電子負荷傳輸方程。

孔相中的氧氣和水蒸氣傳輸方程為：

針對離子聚合物內的水傳輸，膜含水量方程為：

離子聚合物內與孔隙中的水濃度可用Springer關系式[18]描述：

= 0.043 + 17.18? 39.852+ 36.030﹤≤ 1 (7)

其中：d為電滲系數(shù)；e為質子電流；balance為水蒸氣和含水量方程在相界面處因熱力學平衡過程所導致的源項。

電荷方程分為質子勢e方程（13）與電子勢s方程（14）：

模型的主控方程與源項歸納總結于表2中。

在陰極的電化學反應動力學可用Tafel動力學表示：

模型參數(shù)大都來源于MAO等的報道[2]，部分參數(shù)來源于CRC手冊[19]或合作企業(yè)提供的實測結果，材料的物性及部分主要的模型參數(shù)列于表3。

表2 模型主控方程及源項

表3 材料的物性及部分計算參數(shù)

Notes: * data from collaborators; **results from the reconstructed structure.

2.3 初邊界條件和求解方法

當= 0 s時，電極內氧氣濃度為10 mol/m3，孔相及固相中水濃度依據(jù)環(huán)境條件（0.1 MPa，啟動溫度）下相對濕度RH = 38%給定。

氧氣由氣流通道向CL擴散，水蒸氣由CL向氣流通道擴散；在PEM/CL邊界考慮零流條件，則有：

在PEM/CL的邊界上，引入凈水傳輸系數(shù)[20]來衡量催化層與質子交換膜之間水的傳輸量；在CL/GDL界面，膜中水含量傳輸流為0，則有：

其中，/2代表電化學反應產(chǎn)生水的速率。

在= 0和=L處加載電流，=L處的電子電勢固定為0 V，故電荷方程的邊界條件為：

利用商業(yè)計算流體軟件Fluent求解，通過用戶自定義函數(shù)（user-defined functions, UDFs）定義方程的各種源項和物性參數(shù)等；基于有限體積法離散求解，瞬態(tài)項離散采用完全隱式格式，空間導數(shù)項采用中心差分格式差分，采用代數(shù)多重網(wǎng)格（algebraic multi-grid, AMG）迭代法加速收斂。模型中數(shù)值單元為重建結構（圖2）的微元。

3 結果與分析

所求解的算例及其條件如表4所示，電池在恒電流下啟動，首先以基礎算例分析和討論模型結果，然后分析討論變參數(shù)（算例1 ~ 算例4）結果。

表4 算例匯總

3.1 結冰形貌（基礎算例）

模型能預測PEMFC冷啟動過程中冰在催化層內的形貌和分布狀態(tài)。在5 s時開始在孔隙中有冰生成，圖6選擇繪制了不同時刻孔隙中含冰量的三維分布。9.5 s時冰在催化層內已經(jīng)初步形成，隨著時間的推移，冰的份額增加、逐漸覆蓋催化層內孔隙，而且冰占據(jù)的區(qū)域向PEM側擴展。

圖7中截面上的含冰量分布顯示出了冰在孔隙內的局部非均勻生成狀態(tài)。較多冰在孔隙中靠近MS相處生成，孔隙中部位置冰的份額明顯較少，這與結晶過程中晶體會優(yōu)先在固體表面形成的成核理論描述一致。

圖6 催化層內含冰量分布及其隨時間演化

圖7 t =18.5 s時催化層中部x截面的含冰量分布

圖8顯示了不同時刻沿著軸方向的平均含冰量分布，由于重構催化層局部孔隙分布不均勻（圖3）及孔隙內冰的不均勻產(chǎn)生（圖7），冰在占據(jù)孔隙時呈現(xiàn)出不均勻態(tài)勢，而且不均勻程度隨著含冰量的增加愈發(fā)明顯。另外，從圖中可以看出，當含冰量增加時，整體冰含量上升，而冰集中在0.75（=/L，= 1為CL/PEM界面）區(qū)域，并在= 0.5附近出現(xiàn)最大值，而在> 0.75區(qū)間基本沒有冰的生成。冰在催化層內的分布與在= 0和= 1的水傳輸邊界條件及電池的啟動參數(shù)等有關。

圖8 沿x軸方向的含冰量分布及其隨時間演化

3.2 結冰對傳質和電化學反應的影響（基礎算例）

圖9a所示為不同時刻催化層內氧氣平均濃度沿著軸方向的變化情況。從GDL傳輸進來的氧氣由于反應的消耗而沿方向有明顯的下降趨勢，在結冰之前（< 5 s），氧氣在催化層內部的分布趨于穩(wěn)定；當冰開始形成時，由于冰會改變催化層形貌且增加傳質阻力，隨著內部含冰量的增加，氧氣濃度沿方向的下降越來越顯著。盡管如此，可能是由于局部含冰量不均勻分布（參見圖7和圖8），氧氣在催化層內部的傳輸通道并沒有完全被阻斷，在靠近PEM處（= 1）的氧氣濃度沒有因為大量結冰（= 18.5 s）而下降至零。

圖9b所示為離子聚合物內的水含量沿著軸方向的分布。在冷啟動初期（0.5 s與3.5 s），水在催化層內積累，整體含水量上升，當局部達到飽和（= 14）時，水平衡作用使得更多的水析出并于孔隙中發(fā)生相變，而離子聚合物內的水含量保持在飽和值；隨著時間的推移，形成的冰覆蓋界面表面并阻礙了離子聚合物內的水向孔隙內傳輸，導致水含量達飽和值的區(qū)域向PEM側擴展，形成如圖6所示冰的生成區(qū)域向靠近PEM側擴展的現(xiàn)象。

為比較不同時刻催化層內電化學反應的強度以及含冰量對反應的影響，圖10a顯示了6.5 s與18.5 s時的輸運電流密度（c）分布。6.5 s與18.5 s時電化學反應均集中在靠近PEM側（= 1）的區(qū)域，而在陰極CL和陰極GDL的交界面處（= 0），c接近于0；與6.5 s相比，18.5 s時含冰量更高，伴隨著冰覆蓋活性反應面積，電化學反應區(qū)域有向= 0端轉移的趨勢。

結合圖9b與圖10a可以看出，水的分布與電化學反應的主要場所不一致。圖10b給出了6.5 s與18.5 s時的質子電流密度（e）的分布，質子電流密度在= 1附近較大，并且隨著含冰量的增加，質子電流更集中在= 1端附近區(qū)域，可以推斷陰極催化層內電滲作用對水分布的影響不明顯，而圖9b顯示離子聚合物內在= 1附近水含量較低，圖8顯示= 1附近區(qū)域沒有冰生成，可能是由于PEM/CL界面處設置的負的凈水傳輸系數(shù)（見表4）導致的。

圖11 平均過電勢隨含冰量的演化

3.3 電化學有效反應面積與結冰量的關系式

當冰在催化層內形成時，部分活性反應面積將會被覆蓋，影響電化學反應的進行。為了研究電化學有效反應面積與含冰量的關系，定義電化學有效反應面積比eff為催化層內沒有被冰覆蓋的活性反應面積與總活性面積的比值：

圖12示出了各算例下仿真得到的Seff與總含冰量s的關系，并與宏觀冷啟動模型中常采取的線性關系式Seff = 1? s對比，圖中算例1 ~算例4對應表4中算例。

從圖12中可知，仿真結果中的eff均隨著含冰量的增加而減少，各eff曲線的重合度很高（不同算例之間的偏差小于10%），但與常用的線性關系式有明顯偏離；仿真結果的有效反應面積的減少速率更快，覆蓋所有的活性反應面積不需要完全堵塞孔隙。根據(jù)仿真結果擬合得到：

式（24）具有一定的介觀尺度基礎，預期可改進宏觀冷啟動模型的計算精度。

4 結 論

基于PEMFC催化層重建結構建立了電池冷啟動時多場耦合過程的介孔尺度數(shù)值模型。模型考慮了冷啟動時電池內的反應物/產(chǎn)物和電荷傳輸、電化學反應過程，能夠用以深入研究電池冷啟動時陰極催化層內的電?氣?汽（水）傳輸、電化學反應和結冰等。

數(shù)值仿真結果揭示了冷啟動時冰的生成對催化層內電化學反應、水的生成和傳輸、氧氣的傳輸?shù)鹊挠绊懀佻F(xiàn)了冷啟動過程中電極內生成冰的形貌和分布，分析并推演得到電化學活性面積與冰含量的宏觀關系式。該關系式與常用的線性關系式有較大差異，如果用于宏觀模型，預期可以提高模型計算精度。

符號表：

水活度

面積，m2

Bruggeman指數(shù)

摩爾濃度，mol/m3

擴散系數(shù)，m2/s

EW 離子聚合物的當量質量，kg/mol

法拉第常數(shù)，96 485 C/mol

電流負載，A/m2

0交換電流密度，A/m2

e質子電流密度，A/m2

輸運電流密度，A/m3

摩爾質量，kg/mol

反應電子遷移數(shù)

d電滲系數(shù)

速率，s?1

通用氣體常數(shù)，8.314 J/(mol·K)

源項，反應面積比

整體含冰量

溫度，K

0陰極電化學反應熱力學平衡電勢，V

傳輸系數(shù)，凈水傳輸系數(shù)

體積分數(shù)

過電勢，V

離子傳導率，S/m

離子聚合物內含水量

密度，kg/m3

電勢，V

電子傳導率，S/m

角標：

a 環(huán)境

c 陰極

e 電解質

eff 有效的

freeze 結冰

g 氣體

組分

inter 界面

ice 冰

L 總長度

m 模型，離子聚合物

r 反應

ref 參考值

s 固體

sat 飽和值

[1] WEE J H. Applications of proton exchange membrane fuel cell systems[J]. Renewable and sustainable energy reviews,2007, 11(8): 1720-1738. DOI: 10.1016/j.rser.2006.01.005.

[2] MAO L, WANG C Y, TABUCHI Y. A multiphase model for cold start of polymer electrolyte fuel cells[J]. Journal of the electrochemical society, 2007, 154(3): B341-B351. DOI: 10.1149/1.2430651.

[3] JIANG F M, FANG W F, WANG C Y. Non-isothermal cold start of polymer electrolyte fuel cells[J]. Electrochimica acta, 2007, 53(2): 610-621. DOI: 10.1016/j.electacta.2007.07.032.

[4] JIAO K, LI X G. Cold start analysis of polymer electrolyte membrane fuel cells[J]. International journal of hydrogen energy, 2010, 35(10): 5077-5094. DOI: 10.1016/j.ijhydene.2009.09.004.

[5] YAO L, PENG J, ZHANG J B, et al. Numerical investigation of cold-start behavior of polymer electrolyte fuel cells in the presence of super-cooled water[J]. International journal of hydrogen energy, 2018, 43(32): 15505-15520. DOI: 10.1016/j.ijhydene.2018.06.112.

[6] KIM S H, PITSCH H. Reconstruction and effective transport properties of the catalyst layer in PEM fuel cells[J]. Journal of the electrochemical society, 2009, 156(6): B673-B681. DOI: 10.1149/1.3106136.

[7] LANGE K J, SUI P C, DJILALI N. Pore scale modeling of a proton exchange membrane fuel cell catalyst layer: effects of water vapor and temperature[J]. Journal of powersources, 2011, 196(6): 3195-3203. DOI: 10.1016/j.jpowsour.2010.11.118.

[8] WU W, JIANG F M. Microstructure reconstruction and characterization of PEMFC electrodes[J]. International journal of hydrogen energy, 2014, 39(28): 15894-15906. DOI: 10.1016/j.ijhydene.2014.03.074.

[9] WANG G Q, MUKHERJEE P P, WANG C Y. Direct Numerical Simulation (DNS) modeling of PEFC electrodes:part I. Regular microstructure[J]. Electrochimica acta, 2006,51(15): 3139-3150. DOI: 10.1016/j.electacta.2005.09.002.

[10] WANG G Q, MUKHERJEE P P, WANG C Y. Direct Numerical Simulation (DNS) modeling of PEFC electrodes: part II. Random microstructure[J]. Electrochimica acta, 2006, 51(15): 3151-3160. DOI: 10.1016/j.electacta.2005.09.003.

[11] SHEARING P R, BRADLEY R S, GELB J, et al. Using synchrotron X-ray nano-CT to characterize SOFC electrode microstructures in three-dimensions at operating temperature[J]. Electrochemical and solid-state letters, 2011, 14(10): B117-B120. DOI: 10.1149/1.3615824.

[12] SIDDIQUE N A, LIU F Q. Process based reconstruction and simulation of a three-dimensional fuel cell catalyst layer[J]. Electrochimica acta, 2010, 55(19): 5357-5366. DOI: 10.1016/j.electacta.2010.04.059.

[13] 吳偉, 蔣方明, 曾建邦. LiCoO2電池正極微結構重構及有效傳輸系數(shù)預測[J]. 物理化學學報, 2013, 29(11): 2361-2370. DOI: 10.3866/PKU.WHXB201309032.

[14] PROMENTILLA M A B, SUGIYAMA T, HITOMI T, et al. Quantification of tortuosity in hardened cement pastes using synchrotron-based X-ray computed micro- tomography[J]. Cement and concrete research, 2009, 39(6): 548-557. DOI: 10.1016/j.cemconres.2009.03.005.

[15] PATEL K K, PAULSEN J M, DESILVESTRO J. Numerical simulation of porous networks in relation to battery electrodes and separators[J]. Journal of power sources, 2003, 122(2): 144-152. DOI: 10.1016/S0378-7753(03)00399-9.

[16] GE S H, WANG C Y. In situ imaging of liquid water and ice formation in an operating PEFC during cold start[J]. Electrochemical and solid-state letters, 2006, 9(11): A499-A503. DOI: 10.1149/1.2337860.

[17] TAJIRI K, TABUCHI Y, WANG C Y. Isothermal cold start of polymer electrolyte fuel cells[J]. Journal of the electrochemical society, 2007, 154(2): B147-B152. DOI: 10.1149/1.2402124.

[18] SPRINGER T E, ZAWODZINSKI T A, GOTTESFELD S. Polymer electrolyte fuel cell model[J]. Journal of the electrochemical society, 1991, 138(8): 2334-2342. DOI: 10.1149/1.2085971.

[19] HAYNES W M. CRC handbook of chemistry and physics[M]. Raton: CRC Press, 2015.

[20] ZAWODZINSKI T A, DAVEY J, VALERIO J, et al. The water content dependence of electro-osmotic drag in proton-conducting polymer electrolytes[J]. Electrochimica acta, 1995, 40(3): 297-302. DOI: 10.1016/0013-4686(94)00277-8.

[21] WANG C Y. Fundamental models for fuel cell engineering[J]. Chemical reviews, 2004, 104(10): 4727-4766. DOI: 10.1021/cr020718s.

A Mesoscopic Pore-Scale Model of Multi-Disciplinary Coupling in PEMFC Cathode Catalyst Layer during Cold Start

LIAO Zi-hao1,2,3,4, WEI Lin1,2,3, SUO Zhen-bang5, SHI Wei-yu5, JIANG Fang-ming1,2,3

(1. Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences, Guangzhou 510640, China; 2. CAS Key Laboratory of Renewable Energy, Guangzhou 510640, China; 3. Guangdong Provincial Key Laboratory of New and Renewable Energy Research and Development, Guangzhou 510640, China; 4. University of Chinese Academy of Sciences, Beijing 100049, China; 5. Shanghai Hydrogen Propulsion Technology Co., Ltd., Shanghai 201800, China)

In order to have a deep understanding of the cold start of proton exchange membrane fuel cell (PEMFC) from sub-zero temperature, a mesoscopic pore-scale model of multi-disciplinary coupling in cathode catalyst layer (CL) was established. The geometry model is a 3-dimensional pore-structure of CL reconstructed by the random grid method (SGM), and the species transport, charge transport, electrochemical reactions and phase change of water were described by mathematical model. The formation and evolution of ice in CL during cold start were simulated and analyzed, the effects of ice formation and morphology on battery performance were discussed, and a correlation of the electrochemical active reaction area (ECA) to the amount of ice formation in CL was derived from the simulation results.

fuel cells; numerical simulation; mesoscale; catalyst layer; cold start

TK91；TM911.48

A

10.3969/j.issn.2095-560X.2020.02.001

2095-560X（2020）02-0081-10

2019-11-11

2019-12-26

上海汽車工業(yè)科技發(fā)展基金會項目（1706）；廣東省科技發(fā)展專項資金（前沿與關鍵技術創(chuàng)新方向—重大科技專項）項目（2017B01012003）；廣東省自然科學基金項目（2016A030313172）；廣東省自然科學基金?重大基礎研究培育項目（2015A030308019）

蔣方明，E-mail：jiangfm@ms.giec.ac.cn

廖梓豪（1994-），男，博士研究生，主要從事質子交換膜燃料電池數(shù)值模擬及相關實驗研究。

蔣方明（1973-），男，博士，研究員，博士生導師，2001年博士畢業(yè)后曾先后就職于德國的IMM公司、葡萄牙的阿維諾大學、美國的賓州州立大學，2010年底開始回國工作。主要從事綠色能源/動力相關熱物理工程的前沿科學技術問題的研究。