PEM電解槽流場的多物理場耦合建模與分析

童靈華,孫雷波,應(yīng)芳義,姚劍琪,葉夏明,謝長君

(1.國網(wǎng)浙江慈溪市供電有限公司,浙江 慈溪 315300;2.國網(wǎng)浙江省電力有限公司寧波供電公司,浙江 寧波 315100;3.武漢理工大學 自動化學院,武漢 430070)

0 引言

氫儲能具有能量密度高、容量大、壽命長等特點,是實現(xiàn)碳達峰、碳中和目標的重要大規(guī)模儲能技術(shù)[1-4]。作為氫儲能技術(shù)的重要設(shè)備,質(zhì)子交換膜水電解槽(proton exchange membrane water electrolyzer,PEMWE)通過聚合物薄膜傳導(dǎo)質(zhì)子,對功率波動反應(yīng)迅速,與風電、光伏具有良好的匹配性[5]。同時,PEMWE電流密度較高,可顯著提高電解效率,發(fā)展前景良好。

PEMWE的內(nèi)部結(jié)構(gòu)和運行機理復(fù)雜,其性能受多種因素影響,許多學者對其數(shù)值建模進行了研究。Gómez等[6]總結(jié)了過去十年的PEMWE數(shù)值建模研究,認為多數(shù)模型都是基于半經(jīng)驗或經(jīng)驗方程構(gòu)建,高度依賴假設(shè)和近似條件,考慮電解槽動態(tài)操作的模型較少。Nafchi等[7]提出一種一維數(shù)值模型的建立方法,比較了工作溫度、陰極壓力、膜厚度、通道寬度和高度、電流密度等不同參數(shù)對PEMWE性能的影響;Toghyani等[8]采用田口法降低PEMWE的輸入電壓;李建林等[9]采用“組件-單體-陣列-系統(tǒng)”的模型層級分析4種數(shù)值模型的優(yōu)缺點和模型的適用性,并提出PEMWE未來數(shù)值建模的重點研究方向;劉承錫等[10]提出一種可并網(wǎng)的大容量電解槽動態(tài)模型,具備為電網(wǎng)提供快速頻率響應(yīng)的能量。

然而,PEMWE經(jīng)驗建模并不能完全模擬PEMWE內(nèi)部運行過程,特別是流體流動、熱量傳遞等難以簡化表述。近年來,基于多物理場耦合的計算流體動力學(computational fluid dynamics,CFD)模型逐漸應(yīng)用于PEMWE建模。Kaya等[11]利用CFD模型研究了不同負極催化劑對PEMWE性能的影響,發(fā)現(xiàn)Pt-Ir催化劑作為負極催化劑時產(chǎn)氫速率更高。何旭等[12]利用CFD模型分析了多孔運輸層表面接觸面積、進水流速、氧氣生成速率及孔徑大小對氧氣和氫氣泡的生成速率的影響。林楠等[13]建立了一種一維多物理場耦合PEMWE模型,驗證了其電化學特性和傳質(zhì)傳熱特性。

在PEMWE三維CFD建模方面,Zhang等[14]建立了一種單通道PEMWE三維CFD模型研究水流向和流道寬度對電解槽各項指標的影響,認為減小通道深度和增加通道寬度可以改善水傳熱性能。Khatib等[15]建立了3種不同幾何形狀的電解槽三維CFD模型,探究在質(zhì)子交換膜中使用開孔細胞泡沫材料對電解槽性能的影響。一些學者對異形流場的三維CFD建模開展相關(guān)研究,Ojong等[16]創(chuàng)建了一個具有三維多孔傳輸層的半經(jīng)驗全耦合模型,模擬電解槽單電池在多孔金屬棉流場中的反應(yīng)過程。Park等[17]對比通道流場和多孔運輸層(PTL)流場的性能,發(fā)現(xiàn)通道流場性能優(yōu)于PTL流場。Olesen等[18]使用ANSYS對圓平面高壓PEMWE陽極流場中的兩相流進行多物理場耦合分析,發(fā)現(xiàn)減少通道彎曲可以提高電解槽的性能。

在當前PEMWE三維建模研究中,以局部組件建模和單一流場建模為主,缺乏多種流場結(jié)構(gòu)的對比分析,尚未明確不同流場結(jié)構(gòu)與PEMWE電化學性能之間的關(guān)系,因此有必要通過建立PEMWE的三維CFD模型來研究電解槽運行時內(nèi)部單電池的傳熱傳質(zhì)情況,為優(yōu)化流道設(shè)計提供理論依據(jù)。綜上,鑒于目前研究現(xiàn)狀與不足,研究了傳熱、傳質(zhì)和電化學動力學等多物理場耦合下的電解槽工作特性,建立了PEMWE單電池三維CFD模型,仿真分析了平行流場、單通道蛇形流場、雙通道蛇形流場、葉柵形流場等不同流場結(jié)構(gòu)對PEMWE性能的影響。根據(jù)物質(zhì)流速和分布情況,以及電解槽的極化曲線確定PEMWE單電池的最佳流場結(jié)構(gòu)。最后,通過改變溫度、壓力等條件,驗證溫度與壓力對電解槽性能的影響,發(fā)現(xiàn)溫度的升高或壓強的降低都能有效提升電解槽的性能。

1 模型建立

PEMWE單電池由雙極板(bipolar plates,BP)、氣體擴散層(gas diffusion layers,GDL)、催化劑層(catalyst layers,CL)和質(zhì)子交換膜(polymer exchange membrane,PEM)組成。圖1為PEMWE單電池結(jié)構(gòu)示意圖。在陽極,直流電壓將水解離成氧氣、質(zhì)子和電子。質(zhì)子通過質(zhì)子交換膜遷移,與從陰極側(cè)的外部電路傳導(dǎo)的電子結(jié)合形成氫氣。

圖1中的陽極反應(yīng)方程為

2H2O(l)→O2(g)+4H+(aq)+4e-

(1)

陰極反應(yīng)方程為

4H+(aq)+4e-→2H2(g)

(2)

1.1 多物理場模型

在分析PEMWE的多相流情況時,主要考慮電化學動力學、流體動力學和熱傳遞過程。

1.1.1電化學動力學

PEMWE單電池的工作電壓由開路電壓Voc、活化過電壓Vact、濃度差過電壓Vcon和歐姆過電壓Vohm組成,由以下公式表示[19]:

Vcell=Voc+Vact+Vcon+Vohm

(3)

開路電壓由Nernst方程得到:

(4)

式中:R為理想氣體常數(shù);F為法拉第常數(shù);T為工作溫度;pH2和pO 2分別為氫氣和氧氣的分壓;aH2O為電極與膜之間的水活度;E0為電解槽的標準電動勢,E0值受外界環(huán)境的溫度與壓力影響,本文中計算式為

E0=1.229-0.9×10-3(T-298)

(5)

活化極化是由電極反應(yīng)動力學的緩慢性造成的,使用電極表面反應(yīng)的Butler-Volmer方程,可以將活化過電壓用電流密度表示為

(6)

式中:i0表示交換電流密度;αan和αca表示陽極和陰極的電荷轉(zhuǎn)移系數(shù)。

活化過電壓表達式為

(7)

由于電極表面的反應(yīng)物被消耗,電極表面附近的反應(yīng)會與溶液本身產(chǎn)生濃度差,造成電解槽電壓的損失,即濃差極化,故濃度差過電壓的本質(zhì)是粒子擴散的外部表現(xiàn)形式。結(jié)合Nernst方程和Fick定律,濃度差過電壓的表達式為:

(8)

式中:CO2和CH2分別表示膜和電極交界處的氧氣和氫氣濃度;CO2,O和CH2,O分別表示其對應(yīng)參考值。

歐姆過電壓由離子阻力和電極與電極板本身存在的等效電阻引起。根據(jù)歐姆定律,歐姆過電壓與電解過程中的電解電流成線性比例關(guān)系:

Vohm=(2RBP+2RGDL+2RCL+RPEM+Rin)iA

(9)

式中:A表示膜面積。PB和GDL的歐姆電阻是根據(jù)材料的電阻率計算的,而CL和PEM的厚度很薄,其歐姆電阻一般忽略不計。Rin表示離子穿過膜時的電阻,本文中參考文獻[20]相關(guān)研究計算:

(10)

式中:δmem、σmem和λ分別表示質(zhì)子交換膜的厚度、導(dǎo)電率和含水量。

1.1.2流體動力學

質(zhì)量守恒和動量守恒方程用于描述混合物在PEMWE單電池內(nèi)的傳質(zhì)過程[21]:

(11)

式中:ε表示多孔介質(zhì)的孔隙率;μ和ρf表示混合物的黏度和平均密度。Sv表示源項,V表示混合物的體積平均速度[22],其表達式為

(12)

式中:αk和Vk分別為物質(zhì)k的體積分數(shù)和速度。

多孔介質(zhì)中氣液流動的各分量的對流和擴散由Maxwell-Stefan方程表示:

▽·(εVkCk)=▽·(ε1.5Dk▽Ck)+Sk

(13)

式中:Sk是物質(zhì)k的質(zhì)量源項;Ck和Dk分別是物質(zhì)k的摩爾濃度和第k項有效擴散系數(shù),而Dk是關(guān)于溫度和壓力的函數(shù)。

(14)

(15)

氧、氫和水的擴散系數(shù)定義如下:

(16)

1.1.3熱傳遞

PEMWE的熱量來源主要有入口進水帶來的熱量和電流做功時產(chǎn)生的熱量,熱量損失的主要方式為對流散熱、輻射散熱、出口產(chǎn)物以及未反應(yīng)水帶走的熱量。電解槽的熱平衡可以用其能量方程來描述[23]:

(17)

式中:ρeff、Cp,eff、keff表示有效密度、有效熱容和有效熱導(dǎo)率,表達式為

(18)

式中:ρs、Cp,s和ks分別表示固體區(qū)域的密度、熱容和導(dǎo)熱率;ρf、Cp,f和kf分別表示流體混合物的密度、熱容和導(dǎo)熱率。ST是可以通過計算的源項,主要包括不可逆熱Qirrev、電化學反應(yīng)產(chǎn)生的熵熱源Qs、雙極板多孔介質(zhì)和質(zhì)子交換膜產(chǎn)生的歐姆熱源Qohm,表達式為:

(19)

多物理場耦合模型中的參數(shù)見表1。

表1 理論模型的參數(shù)

1.2 三維集合模型

以平行流場為例,圖2為本文中應(yīng)用的PEMWE的幾何模型示意圖,其結(jié)構(gòu)由上到下依次為陽極流場、陽極電極(包括氣體交換層和催化劑)、質(zhì)子交換膜、陰極電極。

圖2 PEMWE單電池幾何模型示意圖

如圖3所示,對3種常見流場形式搭建物理模型,分別是平行流場、單通道蛇形流場和雙通道蛇形流場。另外,設(shè)計葉柵形流場,增大流道與氣體擴散層的接觸面積,探討其電化學性能。平行流場是市面上普遍使用的流場形狀;蛇形流場相較于平行流場更利于反應(yīng)物的接觸,且相對于葉柵形流場具有較少的轉(zhuǎn)角;葉柵形流場具有最大的接觸面積。除流場形態(tài)外,其他幾何參數(shù)相同,物理參數(shù)如表2所示。電極邊長是2 cm,流道的寬度和高度是1 mm,流道間隔是1 mm。

表2 物理參數(shù)

圖3 不同流場形式示意圖

1.3 模型假設(shè)與參數(shù)

為使實驗結(jié)果更直觀,在不影響電解水反應(yīng)過程的前提下,假設(shè):

1) 反應(yīng)物和產(chǎn)物都被認為是不可壓縮的液體或理想氣體,液態(tài)水不會蒸發(fā);

2) 只考慮陽極側(cè)的純水和氧氣混合流體流動,只在陽極側(cè)供水,陰極側(cè)的氫氣產(chǎn)生后直接排出;

3) 電極的多孔介質(zhì)存在各向同性和均勻性;

4) 反應(yīng)物和產(chǎn)物在PEM中的擴散可以忽略不計;

5) 假設(shè)膜含有足夠的水,且電導(dǎo)率僅是溫度的函數(shù)。

2 結(jié)果與討論

2.1 模型驗證

為驗證模型的準確性,將本文模型通過COMSOL Multiphysics仿真計算得到的數(shù)據(jù)與Weiβ等[24]的研究數(shù)據(jù)進行比較,結(jié)果如圖4所示。驗證實驗選取平行流場PEMWE仿真曲線,溫度為25 ℃,壓力為0.1 MPa?？梢钥吹?實驗數(shù)據(jù)與仿真曲線基本擬合,誤差小于3%,證明了模型的有效性。以該模型為基礎(chǔ),討論不同流場對PEMWE性能的影響。

圖4 PEMWE單電池極化曲線(T=25 ℃,P=0.1 MPa)

2.2 流場仿真驗證

2.2.1流速分布

圖5展示了在電流密度為2 A/cm2的情況下,不同流場中電解槽的流速分布情況?？梢钥闯?流場的入口處流體流速最大,隨后流速沿著通道減小。

圖5 不同流場中的速度分布

對比各流場的速度分布,發(fā)現(xiàn)平行流場、單通道蛇形流場、多通道蛇形流場和葉柵形流場中的速度最大值分別為3.94、4.56、4.11和4.34 m/s。平行流場與葉柵形流場的速度峰值主要在進出口,在流道內(nèi)的速度較小,有利于水在陽極充分反應(yīng),且進出口較快的流速易于水的流入和氧氣的流出。單通道蛇形流場在經(jīng)過多個轉(zhuǎn)角后流速有明顯的下降,且流速大于其他2種流場,這會導(dǎo)致反應(yīng)物通過多孔介質(zhì)的遷移較快,雖然較高的流速可以使生成物較快排出,但也會出現(xiàn)反應(yīng)不充分的情況。而雙通道蛇形流場流速分布較為均勻,但在轉(zhuǎn)角處流速從1.5 m/s減小為0.5 m/s,減速情況更為明顯。

2.2.2氧氣濃度分布

在PEM水電解反應(yīng)中,氫氣和氧氣的生成速率主要由電流密度決定,根據(jù)法拉第定律可得:

(20)

式中:n為物質(zhì)的生成速率,由式(20)可知陽極氧氣的生成速率與電流密度成正比。氧氣濃度分布是析氧速率與混合物流動共同作用的結(jié)果。不同流場形態(tài)下,更高的析氧速率表明了更好的電化學性能。同時,流場形態(tài)會影響氧氣在陽極的囤積情況,較高的氧氣濃度會抑制反應(yīng)的持續(xù)進行。

圖6展示了電流密度為2 A/cm2時不同流場電解槽的陽極氧氣濃度分布,各流場從進口端到出口端都有氧氣濃度增大的趨勢。其中,平行流場水電解池的氧氣濃度最大達到了0.63 mol/m3,位于臨近出口一側(cè)的中部,說明水進入多孔介質(zhì)進行了充分反應(yīng),整體上有最高的氧氣生成率,但平行流場存在由于流速較慢、生成物囤積而使氧氣濃度較高的情況;2種蛇形流場的氧氣濃度最大達到0.42和0.41 mol/m3,且從進口到出口均勻增大。同時可以看出,葉柵形流場的氧氣濃度最大達到0.43 mol/m3,主要位于出口一側(cè)的末端和電極的邊緣位置,但并未出現(xiàn)平行流場發(fā)生的氧氣堆積在流道某一部分的情況。

圖6 不同流場中的陽極氧氣分布濃度

2.2.3壓強分布

根據(jù)Navier-Stokes方程中對流體流速與壓強關(guān)系的描述可知,流速越大則壓強越大。結(jié)合前文中對流速的分析,對壓強分布進行分析。

圖7展示了在電流密度為2 A/cm2的情況下不同流場中電解槽的陽極壓力,各流場從進口端到出口端都有壓力減小的趨勢。

圖7 不同流場中的陽極壓力分布

平行流場的壓力減小趨勢幾乎垂直于進出口流道,且減小幅度較為均勻,相鄰流道間的壓力差值約為25 kPa。單通道蛇形流場的壓力減小趨勢平行于長流道,由于單通道蛇形的流速最大,故在此流場形狀下的陽極壓力整體大于其他3種流道形狀,入口處壓力接近300 kPa。雙通道蛇形流場的壓力近似于沿流道方向減小,且每經(jīng)過一個轉(zhuǎn)角壓力降幅約為60 kPa。葉柵形流場的壓力沿對角線減小,此流場由于反應(yīng)較為充分,綜合氣相和液相的壓力減小幅度最大。

2.2.4極化曲線

平行、單通道蛇形、多通道蛇形、葉柵形流場水電解池的極化曲線如圖8所示?？梢钥闯?平行流場中的反應(yīng)物受到流動的限制,導(dǎo)致反應(yīng)點集中,在電流密度為2 A/cm2時極化電壓達1.78 V;2種蛇形流場的性能相差不大,這是因為在中低電流密度時擴散過電位在極化曲線中的占比較小,在電流密度為2 A/cm2時雙通道蛇形流場極化電壓僅比單通道蛇形流場低0.005 V,這是由于雙通道蛇形流場的結(jié)構(gòu)更有利于水在極板通道中流動,不會造成生成物的囤積而導(dǎo)致性能降低;而葉柵形流場中的反應(yīng)物流動性相對最佳,反應(yīng)最充分,在電流密度為2 A/cm2的條件下極化電壓僅有1.72 V。在其他工作條件相同的情況下,制取相同量的氫氣,葉柵形流場的PEMWE消耗的電能最少,即電解槽性能最佳。

圖8 不同流場下的極化曲線

2.3 參數(shù)敏感性分析

對于同一PEMWE單電池結(jié)構(gòu),工作溫度和堆內(nèi)壓力是影響其電化學性能的主要參數(shù)。選取20、40、60和80 ℃這4個不同溫度分別計算壓強為0.4 MPa條件下的極化曲線,如圖9所示?？梢钥闯?當PEMWE工作溫度較高時,極化電壓更低。這是由于高溫能促進反應(yīng)離子的擴散和向陰極的遷移,因此較高的溫度可以增強質(zhì)子交換膜的傳輸特性,降低電解槽所需電壓。在電流密度為2 A/cm2時,其極化電壓分別為1.727、1.721、1.714和1.707 V。

圖9 不同溫度下的極化曲線

選取0.1、0.2、0.4和0.8 MPa這4個不同壓強在80 ℃工作溫度下的極化曲線,如圖10所示?？梢钥闯?陰極室的壓力增大,阻礙氫離子向陰極移動,從而阻礙電解反應(yīng)的進行,因此電解電壓會隨著壓力的升高而升高。在電流密度為2 A/cm2時,其極化電壓分別為1.686、1.696、1.707和1.717 V。

圖10 不同壓強下的極化曲線

3 結(jié)論

1) PEMWE電氫轉(zhuǎn)化率受到流場結(jié)構(gòu)的影響,其中平行流場流速為3.94 m/s,相較于其他流場速度較慢,析氧反應(yīng)在膜電極中分布不均勻,產(chǎn)生了較大的電化學與濃差過電位,在電流密度為2 A/cm2時極化電壓達到1.78 V。在單通道蛇形流場和雙通道蛇形流場中流速較快,但過快的流速導(dǎo)致反應(yīng)不充分,且壓力減小較為緩慢。葉柵形流場出口流速較快,為4.34 m/s,反應(yīng)充分,產(chǎn)生的電動勢低于平行流場3.5%、低于蛇形流場2.3%,壓力減小幅度最大,性能相對最好。

2) 在葉柵形流場的條件下,溫度每升高20 ℃,電解槽極化電壓下降0.4%;槽內(nèi)壓力每增大1倍,電解槽極化電壓上升0.59%,即升高電解槽溫度或降低反應(yīng)室壓力都有助于電解反應(yīng)的進行,提升PEMWE的性能。