PEM電解槽流場(chǎng)的多物理場(chǎng)耦合建模與分析

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

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

0 引言

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

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

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

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

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

1 模型建立

PEMWE單電池由雙極板(bipolar plates,BP)、氣體擴(kuò)散層(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 多物理場(chǎng)模型

在分析PEMWE的多相流情況時(shí),主要考慮電化學(xué)動(dòng)力學(xué)、流體動(dòng)力學(xué)和熱傳遞過程。

1.1.1電化學(xué)動(dòng)力學(xué)

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

Vcell=Voc+Vact+Vcon+Vohm

(3)

開路電壓由Nernst方程得到:

(4)

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

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

(5)

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

(6)

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

活化過電壓表達(dá)式為

(7)

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

(8)

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

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

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

(9)

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

(10)

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

1.1.2流體動(dòng)力學(xué)

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

(11)

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

(12)

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

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

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

(13)

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

(14)

(15)

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

(16)

1.1.3熱傳遞

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

(17)

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

(18)

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

(19)

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

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

1.2 三維集合模型

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

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

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

表2 物理參數(shù)

圖3 不同流場(chǎng)形式示意圖

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

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

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

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

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

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

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

2 結(jié)果與討論

2.1 模型驗(yàn)證

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

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

2.2 流場(chǎng)仿真驗(yàn)證

2.2.1流速分布

圖5展示了在電流密度為2 A/cm2的情況下,不同流場(chǎng)中電解槽的流速分布情況。可以看出,流場(chǎng)的入口處流體流速最大,隨后流速沿著通道減小。

圖5 不同流場(chǎng)中的速度分布

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

2.2.2氧氣濃度分布

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

(20)

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

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

圖6 不同流場(chǎng)中的陽極氧氣分布濃度

2.2.3壓強(qiáng)分布

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

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

圖7 不同流場(chǎng)中的陽極壓力分布

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

2.2.4極化曲線

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

圖8 不同流場(chǎng)下的極化曲線

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

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

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

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

圖10 不同壓強(qiáng)下的極化曲線

3 結(jié)論

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

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