摘要：針對(duì)質(zhì)子交換膜燃料電池在機(jī)械應(yīng)力下的氣-液兩相流進(jìn)行數(shù)學(xué)模擬研究，建立了一個(gè)二維質(zhì)子交換膜燃料電池非等溫兩相流多物理場(chǎng)穩(wěn)態(tài)模型.該模型綜合考慮了固體力學(xué)、電化學(xué)、傳熱傳質(zhì)以及氣液兩相流的物理因素，研究了質(zhì)子交換膜燃料電池在機(jī)械應(yīng)力作用下的兩相流分布.計(jì)算結(jié)果顯示：在機(jī)械應(yīng)力作用下，燃料電池肋板下方的多孔介質(zhì)應(yīng)力明顯大于流道下方的應(yīng)力，且在肋板和流道交界處下方的氣體擴(kuò)散層會(huì)產(chǎn)生明顯的應(yīng)力集中現(xiàn)象;隨著電流密度的增加，陰極相對(duì)濕度逐漸增加，但陽(yáng)極相對(duì)濕度會(huì)減小;液態(tài)水僅在陰極產(chǎn)生且主要在肋板下方的多孔介質(zhì)內(nèi)形成，其在陰極的飽和度隨電流密度的增加而不斷增加.

關(guān)鍵詞：質(zhì)子交換膜燃料電池;機(jī)械應(yīng)力;多物理場(chǎng);兩相流;液態(tài)水

中圖分類號(hào)：TM911.4文獻(xiàn)標(biāo)志碼：A

Working Parameters Distribution of Fuel Cell under Mechanical Stress

ZHANG Heng，CHEN Ben，SUI Bangjie

（School of Automotive Engineering，Wuhan University of Technology，Wuhan 430070，China）

Abstract：The mathematical simulation of gas-liquid two-phase flow of the proton exchange membrane （PEM）fuel cell under mechanical stress was studied. In this study，a two-dimensional，non-isothermal two-phase flow Multiphysics steady-state model of PEM fuel cells was established. The model comprehensively considered the solid mechanics，electrochemistry，heat and mass transfer and gas-liquid two-phase flow. The two-phase flow distribution of PEM fuel cells under mechanical stress was studied. The computational results showed that the stress of the porous medium under ribs was significantly greater than the stress under the flow channel. Stress concentration obviously occurred at the junction of the ribs and the flow channel. Liquid water was only condensed at the cathode and mainly formed in the porous medium under the ribs. As the current density increased，the cathode relative humidity gradually increased，however，the relative humidity of the anode decreased. Cathode liquid water saturation increased when current density increased.

Key words：proton exchange membrane fuel cell;mechanical stress;Multiphysics;two-phase;liquid water

質(zhì)子交換膜燃料電池由于具有工作效率高、能量密度高和污染低等優(yōu)勢(shì)被視為目前最有前景的可替代清潔能源[1].為提高燃料電池的性能，全面了解電池內(nèi)部的各種傳輸現(xiàn)象是十分必要的.燃料電池包含5個(gè)關(guān)鍵部件，從外到內(nèi)分別為雙極板（Bipolar plate，BPP）、氣體擴(kuò)散層（Gas diffusion layer，GDL）、微孔層（Micro-porous layer，MPL）、催化層（Catalyst layer，CL）和質(zhì)子交換膜（Proton exchange membrane，PEM）[2].由于單電池所產(chǎn)生的電能有限，為了滿足大功率的使用要求，通常將多個(gè)單電池進(jìn)行串聯(lián)形成燃料電池堆，以便輸出足夠的功率.在燃料電池堆的裝配過(guò)程中，過(guò)小的裝配力會(huì)導(dǎo)致反應(yīng)氣體泄漏和產(chǎn)生較大的導(dǎo)熱導(dǎo)電阻力[3];過(guò)大的裝配力會(huì)導(dǎo)致電池內(nèi)部產(chǎn)生較高的氣體傳輸阻力，阻礙反應(yīng)氣體的傳輸和產(chǎn)物水的排出，而且還會(huì)導(dǎo)致機(jī)械損壞[4].在工程實(shí)際中，一般選用使GDL產(chǎn)生20%形變量（載荷約為1.8 MPa）的機(jī)械力來(lái)裝配燃料電池堆，此時(shí)燃料電池堆能達(dá)到最佳的工作性能[5].Millichamp等人[3]描述了機(jī)械壓縮對(duì)燃料電池中組件尺寸的影響，并總結(jié)了實(shí)現(xiàn)燃料電池堆最佳性能的機(jī)械壓力.Chi等人[6]用實(shí)驗(yàn)和理論方法研究夾緊力導(dǎo)致的GDL孔隙率不均勻性對(duì)燃料電池性能的影響. Zhou等人[7]建立了一個(gè)三維多相模型分析不均勻壓縮對(duì)燃料電池內(nèi)部傳輸參數(shù)以及對(duì)電池傳輸過(guò)程的影響.這些模型考慮了機(jī)械應(yīng)力，但是沒(méi)有考慮組件內(nèi)部明顯的機(jī)械形變.由于GDL的楊氏模量較其他組件很低，所以在機(jī)械應(yīng)力下，GDL為最容易產(chǎn)生機(jī)械形變的組件[8].與此同時(shí)，不同的兩相流模型被建立來(lái)研究燃料電池中的液態(tài)水傳輸.Ince等人[9]基于XCT實(shí)驗(yàn)，研究了燃料電池在不同的壓力下的水傳輸特征，結(jié)果表明GDL中液態(tài)水飽和度在壓縮的碳纖維區(qū)域比未壓縮的區(qū)域更高.Havaej等人[10]建立了一個(gè)三維瞬態(tài)的兩相流模型模擬陰極半電池的性能.Caulk等人[11]用兩相流模型研究了極限電流和液態(tài)水在GDL中的傳輸阻力.Vetter等人[12]用一維非等溫兩相流MATLAB模型分析了從流道到催化層的液態(tài)水飽和度分布情況.但這些研究沒(méi)有考慮機(jī)械應(yīng)力對(duì)液態(tài)水分布的影響.

為了完善現(xiàn)有的兩相流計(jì)算模型，本文首先求解單電池在機(jī)械應(yīng)力作用下的應(yīng)力分布，并將GDL 內(nèi)與應(yīng)力形變有關(guān)的傳輸參數(shù)納入新的幾何模型中，得到了實(shí)際的傳輸性能參數(shù)[13-14]（如孔隙滲透率、電導(dǎo)率、熱導(dǎo)率、接觸電阻等）.然后，研究在機(jī)械應(yīng)力作用下不同電流密度的液態(tài)水分布規(guī)律.此外，模型還考慮了液態(tài)水在膜中反擴(kuò)散和電拖曳兩種傳輸方式，從而能全面求解液態(tài)水在各個(gè)組件中的分布.該文的研究對(duì)燃料電池的水管理有非常重要的指導(dǎo)作用.在工程實(shí)踐中，可以對(duì)電池內(nèi)部的水分布做出精確的預(yù)測(cè)，有利于提高工程效率.

1模型建立

本計(jì)算模型分為兩部分，包括圖1所示的未受機(jī)械應(yīng)力時(shí)的計(jì)算域和圖2所示的受到機(jī)械應(yīng)力之后發(fā)生形變的計(jì)算域.兩部分的計(jì)算域組成部件從外及內(nèi)為陰陽(yáng)極BPP、GDL、MPL、CL和PEM，以及BPP內(nèi)部流道：陽(yáng)極流道注入氫氣，陰極流道注入空氣.從圖2可以看出GDL受到機(jī)械應(yīng)力后發(fā)生了明顯的形變且出現(xiàn)了入侵流道的現(xiàn)象[25-16].

1.1模型假設(shè)

本研究中，建立了一個(gè)二維穩(wěn)態(tài)的PEM燃料電池兩相流模型.為了簡(jiǎn)化計(jì)算，并保證計(jì)算的可行性，本模型基于以下假設(shè)條件[13]：

1）本模型計(jì)算是處于穩(wěn)態(tài);

2）反應(yīng)氣體均為理想氣體;

3）MPL、CL以及PEM的應(yīng)力形變被忽略，且傳輸參數(shù)均為各向同性，孔徑分布均勻;

4）GDL里面的形變?yōu)閺椥孕巫?

5）機(jī)械應(yīng)力里面的濕熱效應(yīng)被忽略.

1.2控制方程

GDL在模型計(jì)算中視為可壓縮的彈性材料，機(jī)械應(yīng)力形變方程為胡克定律：

電子傳輸方程為歐姆定律：

式中：σ_e為電子傳導(dǎo)率;?_e為電子電勢(shì)差;j_e為電子通量;i為電流密度.

質(zhì)子傳輸方程為歐姆定律：

式中：σ_P為質(zhì)子傳導(dǎo)率;?_P為質(zhì)子電勢(shì)差;j_P為質(zhì)子通量.

熱傳導(dǎo)方程如下：

式中：K為熱導(dǎo)率;T為溫度;S_T，e為電子流動(dòng)產(chǎn)生的熱源;S_T，p為質(zhì)子流動(dòng)產(chǎn)生的熱源;S_T，r為電化學(xué)反應(yīng)熱源.

氣體傳質(zhì)方程為Maxwell-Stefan方程：

式中：ρ為密度;ω_i，ω_j為質(zhì)量分?jǐn)?shù);D_ij為二元?dú)怏w的擴(kuò)散系數(shù);M為摩爾質(zhì)量;x_i，x_j為摩爾質(zhì)量分?jǐn)?shù);p為壓強(qiáng);S_F為電化學(xué)反應(yīng)源項(xiàng);S_ec為蒸發(fā)/冷凝源項(xiàng);S_ad為吸附/解析源項(xiàng).

液態(tài)水的傳輸方程服從達(dá)西定律[17]：

式中：κ為滲透率;μ為液體水的黏度;V_W為水的摩爾體積;P_C為毛細(xì)壓力;S為液態(tài)水的飽和度.

氣體擴(kuò)散介質(zhì)的毛細(xì)壓力和飽和度的關(guān)系[12]：

式中：C為總間隙氣體濃度;x_sat為飽和摩爾分?jǐn)?shù);P_sat為飽和壓強(qiáng).

1.3形變修正方程

式中[6]：CR為壓縮率;V_comp，為GDL受到機(jī)械應(yīng)力壓縮形變之后的體積;V為氣體擴(kuò)散層受到壓縮之前的初始體積.

式中：ε_comp為受到機(jī)械應(yīng)力之后的孔隙率;ε₀為初始孔隙率.

式中：τ_comp為壓縮之后的曲率.

式中：x_comp為滲透率;其中C_k為Kozeny常數(shù)（C_k=6）;D_r為氣體擴(kuò)散介質(zhì)中的平均孔直徑（D_r=7×10^-6m）.

式中[18]：r為氣體擴(kuò)散層與雙極板之間的接觸電阻率（mΩ·cm^-2）;p為雙極板與擴(kuò)散層之間的接觸壓力（Pa）;A，B，C 均為常數(shù)（A=3.32 mΩ·cm^- B=1.01 MPa，C=0.534）.

1.4液態(tài)水飽和度修正方程[14]

式中：S為液態(tài)水的飽和度;ε₁為考慮液態(tài)水之后的孔隙率;τ₁為考慮液態(tài)水飽和度之后的曲率;K₁為液態(tài)水飽和度修正之后的熱導(dǎo)率;χ₁為液態(tài)水飽和度修正之后的滲透率.

1.5源項(xiàng)

表1為本模型計(jì)算公式的源項(xiàng)匯總，包含多孔介質(zhì)層（GDL、MPL、CL）和PEM.

在電池的陽(yáng)極氫氣發(fā)生氧化反應(yīng)產(chǎn)生質(zhì)子和電子，此為電子和質(zhì)子的源項(xiàng).在電池的陰極氧氣和陽(yáng)極產(chǎn)生的電子和質(zhì)子發(fā)生還原反應(yīng)產(chǎn)生水.在本文中，假定反應(yīng)產(chǎn)生的水均為氣態(tài)水S_F.

氧氣的消耗為S_F/ 氫氣的消耗為S_F.S_ad為在CL里面液態(tài)水和離聚物形態(tài)的水轉(zhuǎn)換的源項(xiàng)：

式中：κ_a和κ_d分別為吸附和解吸質(zhì)量傳輸系數(shù);L為催化層厚度;V_m為離聚物摩爾體積，λ_eq為離聚物中水的等效含量，λ為離聚物的水含量.

離聚物中水的等效含量λ_eq表達(dá)式為[19]：

氣相水和液相水的相變轉(zhuǎn)化，由以下轉(zhuǎn)變方程式?jīng)Q定[12]：

式中：σ_e為電子電導(dǎo)率;σ_p，為質(zhì)子電導(dǎo)率;φ_e為電子電勢(shì);φ_p為質(zhì)子電勢(shì);i為電流密度;E⁰為等效電勢(shì).

1.6邊界條件

表2總結(jié)了各個(gè)物理場(chǎng)的邊界條件.該模型中所有的反應(yīng)氣體，假定均不能滲透過(guò)質(zhì)子交換膜，且質(zhì)子交換膜可以阻止所有的電子通過(guò).在模型中，通過(guò)控制相對(duì)濕度來(lái)控制摩爾分?jǐn)?shù)，其中相對(duì)濕度（RH）和摩爾分?jǐn)?shù)的關(guān)系為：

2結(jié)果和分析

2.1模型驗(yàn)證

圖3所示的極化曲線，將實(shí)驗(yàn)數(shù)據(jù)和模型數(shù)據(jù)進(jìn)行對(duì)比，本實(shí)驗(yàn)數(shù)據(jù)陽(yáng)極和陰極在相同的工作條件下測(cè)得（陽(yáng)極和陰極工作條件一樣：壓強(qiáng)=200 kPa、溫度=353. 15 K、相對(duì)濕度=80%），實(shí)驗(yàn)數(shù)據(jù)和模型仿真數(shù)據(jù)吻合良好，證明了該模型的合理性.本組實(shí)驗(yàn)的特征參數(shù)總結(jié)見(jiàn)表3.

2.2PEM燃料電池的力學(xué)分析

圖4為PEM燃料電池在受到裝配力時(shí)，多孔介質(zhì)的應(yīng)力分布情況（文中圖形結(jié)果僅為多孔介質(zhì)區(qū)域）.從圖4的應(yīng)力分布可以發(fā)現(xiàn)：肋板下方的多孔介質(zhì)受到的應(yīng)力遠(yuǎn)大于流道下方的多孔介質(zhì)，且流道下面的GDL會(huì)產(chǎn)生入侵流道的現(xiàn)象.與此同時(shí)，在肋板和流道交界處下方的GDL有明顯的應(yīng)力集中現(xiàn)象，此處的應(yīng)力值為多孔介質(zhì)中的最大值5 MPa.圖5為燃料電池受到機(jī)械應(yīng)力后的位移分布，從圖中可以看出：由于肋板和BPP直接接觸，位于肋板下方GDL的形變量要大于位于流道下方GDL的形變量，流道下方GDL的形變量幾乎為0. MPL、CL以及PEM的形變量相對(duì)于GDL的形變量而言是極小的.

2.3PEM燃料電池的溫度分布

圖6為PEM燃料電池在不同電流密度下（i=0.17 A·cm^-2、1.05 A·cm^-2、1.66 A·cm^-2）的溫度分布.從圖6可以看出：流道下方多孔介質(zhì)的溫度高于肋板下方多孔介質(zhì)的溫度.由于機(jī)械壓縮，肋板下方多孔介質(zhì)的孔隙率相較于流道下方多孔介質(zhì)的孔隙率要小，這樣導(dǎo)致肋板下方的氣體傳輸阻力要大于流道下方的氣體傳輸阻力，從而使流道下方的多孔介質(zhì)內(nèi)部反應(yīng)氣體的濃度較高，相應(yīng)的電化學(xué)反應(yīng)產(chǎn)生的熱量較多，溫度也較高.還可以發(fā)現(xiàn)陰極的溫度要高于陽(yáng)極，這是由于在PEM燃料電池中，陽(yáng)極和陰極是分別發(fā)生兩個(gè)半電化學(xué)反應(yīng)，陰極的還原反應(yīng)的反應(yīng)焓要大于陽(yáng)極氧化反應(yīng)的反應(yīng)焓，所以陰極的溫度要高于陽(yáng)極的溫度.隨著電流密度的增加，電池內(nèi)部的溫度會(huì)增加，陽(yáng)極和陰極的溫差會(huì)隨著電流密度的增加而增加.

2.4PEM燃料電池的相對(duì)濕度分布

圖7是PEM燃料電池在三個(gè)不同的電流密度下的（i=0.17 A·cm^-2、1.05 A·cm^-2、1.66 A·cm^-2）相對(duì)濕度分布.從圖7（注：陽(yáng)極和陰極的圖例不同）可以發(fā)現(xiàn)：陰極的相對(duì)濕度隨著電流密度的增加而逐漸增加，并且會(huì)超過(guò) 尤其是在肋板下方的多孔介質(zhì)區(qū)域增加明顯.這是因?yàn)槔甙逑路降亩嗫捉橘|(zhì)在機(jī)械應(yīng)力的作用下，傳輸參數(shù)如孔隙率和滲透率的減小，傳輸阻力的增大使得電化學(xué)產(chǎn)生的水不容易排出，加上肋板下方的多孔介質(zhì)區(qū)域溫度要低于流道下方的多孔介質(zhì)區(qū)域溫度，所以肋板下方區(qū)域的相對(duì)濕度要高于流道下方的.陽(yáng)極的相對(duì)濕度分布則與陰極情況相反，陽(yáng)極的相對(duì)濕度隨著電流密度的增加而逐漸減小，且陽(yáng)極的相對(duì)濕度均小于1.這是因?yàn)樗牧魍渴怯呻娡弦泛头磾U(kuò)散兩種因素共同決定的.隨著電流密度的增加，電拖拽作用增強(qiáng)，陽(yáng)極的水會(huì)隨著質(zhì)子向陰極流動(dòng)，當(dāng)電拖拽作用的流通量大于反擴(kuò)散作用的流通量時(shí)，陽(yáng)極內(nèi)的水會(huì)隨著電流密度的增加而減少.此外，陰極流道下方的水含量大于肋板下方的水含量，使得流道下方區(qū)域反擴(kuò)散的作用大于肋板下方區(qū)域的反擴(kuò)散作用，所以陽(yáng)極流道下方的相對(duì)濕度大于肋板下方的相對(duì)濕度.隨著電流密度的增加，CL上的相對(duì)濕度首先在其兩邊變大，然后朝著CL中間區(qū)域增加.當(dāng)電流密度較小時(shí)，CL中間區(qū)域的反應(yīng)物濃度較高，使其中間區(qū)域溫度較高，所以兩邊區(qū)域的相對(duì)濕度高于中間區(qū)域的相對(duì)濕度.當(dāng)電流密度繼續(xù)增加時(shí)，反應(yīng)產(chǎn)生的水占了主導(dǎo)因素.因而，CL中間區(qū)域的相對(duì)濕度隨著水濃度的明顯增加而變大.

2.5PEM燃料電池液態(tài)水飽和度分布

從圖7可以看出：在不同電流密度的情況下，由于陽(yáng)極的相對(duì)濕度均小于 所以陽(yáng)極一直沒(méi)有液態(tài)水產(chǎn)生，故在圖8中沒(méi)有繪制陽(yáng)極液態(tài)水飽和度的分布.圖8是PEM燃料電池在三個(gè)不同的電流密度下的（i=0.17 A·cm^-2、1.05 A·cm^-2、1.66 A·cm^-2）陰極液態(tài)水飽和度分布.陰極的相對(duì)濕度會(huì)在電流密度達(dá)到一定數(shù)值的時(shí)候大于 此時(shí)會(huì)產(chǎn)生液態(tài)水.且陰極液態(tài)水飽和度隨著電流密度的增加而逐漸增加[ 20].當(dāng)電流密度為1.66 A·cm^-2的時(shí)候，液態(tài)水飽和度可以達(dá)到0.3.液態(tài)水飽和度在肋板下方的多孔介質(zhì)區(qū)域要大于流道下方，這是因?yàn)槔甙逑路降目紫堵屎蜐B透率要小于流道下方的，所以產(chǎn)生的液態(tài)水更容易在肋板下方積累.且由于GDL區(qū)域的溫度要遠(yuǎn)小于CL區(qū)域的溫度，液態(tài)水飽和度的最大值會(huì)出現(xiàn)在肋板下方遠(yuǎn)離CL的GDL區(qū)域.

3結(jié)論

本模型結(jié)合了機(jī)械應(yīng)力和傳熱傳質(zhì)、電化學(xué)的多物理場(chǎng)模型求解氣-液兩相流，研究了機(jī)械應(yīng)力下PEM燃料電池的氣-液兩相流工作參數(shù)分布情況.本研究結(jié)論如下：

1）燃料電池在受到裝配機(jī)械應(yīng)力時(shí)，肋板下方的GDL形變量最大，其他部件的形變量相對(duì)于GDL 形變幾乎可以忽略.肋板下方的多孔介質(zhì)受到的機(jī)械應(yīng)力要大于流道下方多孔介質(zhì)受到的應(yīng)力，且在肋板與流道交界處GDL有應(yīng)力集中現(xiàn)象.

2）燃料電池陰極的溫度要高于陽(yáng)極，陰極CL的溫度最高，隨著電池的電流密度增加，電池的溫度將增加，溫差也會(huì)變大.當(dāng)電流密度為1.66 A·cm^-2時(shí)，溫差為10 K.

3）隨著電流密度的增加，陰極的相對(duì)濕度從初始的0.8增加到1. 但陽(yáng)極相對(duì)濕度從初始的0.8減小到0.6左右.

4）陰極隨著電流密度的增加會(huì)逐漸產(chǎn)生液態(tài)水，電流密度越大，液態(tài)水的飽和度越高.當(dāng)電流密度為1.66 A·cm^-2時(shí)，飽和度為0.3.肋板下方多孔介質(zhì)層中的液態(tài)水飽和度要大于流道下方的多孔介質(zhì)層，且飽和度的最大值會(huì)發(fā)生在離肋板最近的GDL 區(qū)域.

參考文獻(xiàn)

[1] MAJLAN E H，ROHENDI D ，DAUD W R W，et al. Electrode for proton exchange membrane fuel cells：a review [J]. Renewable and Sustainable Energy Reviews ，2018，89：117-134 .

[2] PRIYA K，SATHISHKUMAR K，RAJASEKAR N. A comprehensive review on parameter estimation techniques for Proton Exchange Membrane fuel cell modelling[J] . Renewable and Sustainable Energy Reviews，2018，93：121-144 .

[3] MILLICHAMP J，MASON T J，NEVILLE T P，et al. Mechanisms and effects of mechanical compression and dimensional change in polymer electrolyte fuel cells-A review[J] . Journal of Power Sources，2015，284：305-320.

[4] SIMON C，HASCHEF，GASTEIGER H A. Influence ofthe gas diffusion layer compression on the oxygen transport in PEM fuel cells at high water saturation levels[J] . Journal of the Electrochemical Society，2017，164（6）：F591-F599.

[5] DAFALLA A M，JIANG F M. Stresses and their impacts on proton exchange membrane fuel cells：a review[J] . International Journal of Hydrogen Energy ，2018，43（4）：2327-2348.

[6] CHI P H，CHAN S H，WENG F B ，et al. On the effects of nonuniform property distribution due to compression in the gas diffusion layer of a PEMFC[J] . International Journal of Hydrogen Energy，2010，35（7）：2936-2948.

[7] ZHOU Y B，JIAO K，DU Q，et al. Gas diffusion layer deformation and its effect on the transport characteristics and performance of proton exchange membrane fuel cell[J] . International Journal of Hydrogen Energy ，201 38（29）：12891-12903 .

[8] ZHANG H，RAHMAN M A，MOJICA F，et al. A comprehensive two-phase proton exchange membrane fuel cell model coupled with anisotropic properties and mechanical deformation of the gas diffusion layer[J] . Electrochemistry Acta，202 382：138273 .

[9] INCE U U，MARKOTTER H，GEORGE M G，et al. Effects of compression on water distribution in gas diffusion layer materials of PEMFC in a point injection device by means of synchrotron X-ray imaging[J] . International Journal of Hydrogen Energy，2018，43 （1）：391-406.

[10] HAVAEJ P. A numerical investigation of the performance of Polymer Electrolyte Membrane fuel cell with the converging-diverging flow field using two-phase flow modeling[J] . Energy，2019，182：656-672.

[11] CAULK D A，BAKER D R.Modeling two-phase water transport in hydrophobic diffusion media for PEM fuel cells[J] . Journal of the Electrochemical Society，201 158（4）：B384.

[12] VETTER R，SCHUMACHER J O.Free open reference implementation of a two-phase PEM fuel cell model[J] . Computer Physics Communications，2019，234：223-234.

[13] ZHANG H，XIAO L S，CHUANG P Y A，et al. Coupled stress-strain and transport in proton exchange membrane fuel cell with metallic bipolar plates[J] . Applied Energy，2019，251：113316.

[14] WEBER A Z，NEWMAN J.Modeling transport in polymer-electrolyte fuel cells[J] . Chemical Reviews，200 104（10）：4679-4726.

[15] LI S A，SUNDEN B. Effects of gas diffusion layer deformation on the transport phenomena and performance of PEM fuel cells with interdigitated flow fields[J] . International Journal of Hydrogen Energy，2018，43（33）：16279-16292.

[16] WANG J T，YUAN J L，SUNDEN B. On electric resistance effects of non-homogeneous GDL deformation in a PEM fuel cell[J]. International Journal of Hydrogen Energy，2017，42（47）：28537- 28548.

[17] WEBER A Z，BORUP R L，DARLING R M，et al. A critical review of modeling transport phenomena in polymer-electrolyte fuel cells[J]. Journal of the Electrochemical Society，201 161（12）：F1254-F1299.

[18] ZHOU P，WU C W，MA G J. Contact resistance prediction and structure optimization of bipolar plates [J]. Journal of Power Sources，2006，159（2）：1115-1122.

[19] VIJAYARAGHAVAN K，DEVAAL J，NARIMANI M. Dynamic model of oxygen starved proton exchange membrane fuel-cell using hybrid analytical-numerical method[J]. Journal of Power Sources，2015，285：291-302.

[20] MUIRHEAD D，BANERJEE R，GEORGE M G，et al. Liquid water saturation and oxygen transport resistance in polymer electrolyte membrane fuel cell gas diffusion layers [J] . Electrochemistry Acta，2018，274：250-265.