面向燃料電池雙極板的不銹鋼表相摻雜耐蝕性高通量計算與分析

張浩明, 張 頔, 徐竹田, 彭林法,來新民

(1.上海交通大學 上海市復雜薄板結(jié)構(gòu)數(shù)字化制造重點實驗室,上海 200240;2.上海交通大學 機械系統(tǒng)與振動國家重點實驗室,上海 200240)

1 引 言

質(zhì)子交換膜燃料電池(Proton exchange membrane fuel cells，PEMFCs)由于其能量密度高、響應(yīng)快速、環(huán)境友好等特點，被認為是具有廣泛的應(yīng)用前景和極高的研究價值的新一代氫能源轉(zhuǎn)換裝置，并且在交通運輸、航空航天、能源儲備等各方面取的廣泛應(yīng)用[1，2]. 雙極板(Bipolar plate，BPP)是PEMFCs電堆中的關(guān)鍵部件，承擔著支撐反應(yīng)空間、分配反應(yīng)氣體和冷卻水、傳導電流等作用. 近些年超薄、低成本、易于大批量制造的不銹鋼雙極板成為燃料電池產(chǎn)業(yè)化的研究熱點[3，4].

燃料電池極板工作于酸性(pH≈2-5)、高溫(約70-80℃)、含有F-等腐蝕性離子、受電極電勢作用的腐蝕性溶液環(huán)境，眾多研究表明不銹鋼材料受到腐蝕可能會造成接觸電阻上升、離子析出(Fe、Cr離子等)毒害質(zhì)子交換膜等問題[5]，因此對雙極板耐蝕性提出了更高要求. 目前常用的提高不銹鋼材料性能的途徑是采用PVD[6]、CVD[7]、PIRAC[8]、電沉積[9，10]等方法進行涂層處理，但是采用涂層改性會給雙極板的生產(chǎn)帶來額外的成本，而且涂層本身也存在開裂、脫附等失效的可能性. 因此也有學者進行了免涂層不銹鋼基材的研究，通過對基材進行離子注入、成分調(diào)控等改性方式來使其滿足PEMFCs性能需求，例如Feng在316L不銹鋼中采用離子注入的方法注入Ni-Cr或者Ag，從而提高材料耐蝕、導電性能[11，12]；Huang設(shè)計團簇式為[Ni-Fe13-xNix-1]Cr3的不銹鋼體系，通過多組實驗得出x約為3時具有最適合PEMFCs雙極板基材的性能[13].

相較于經(jīng)驗指導、實驗驗證的傳統(tǒng)方法，理論計算與模擬常常可以從更小的尺度上探究材料的機理并指導材料的開發(fā)，進而降低試錯成本、縮短研發(fā)周期[14]. 不少學者已經(jīng)開始使用理論計算方法對新型不銹鋼設(shè)計展開研究，Vitos采用第一性原理方法計算了Fe100-(m+n)CrmNin合金及其加入多種金屬元素后的剪切模量和體積模量，認為主元素配比為Fe58Cr18Ni24時具有適中的硬度的同時具有較好的延展性和耐蝕性，且加入少量Os或者Ir元素能進一步提高這些性能[15]；Feng基于密度泛函理論(DFT)計算和加壓冶金，設(shè)計了高N含量(>0.33wt%)的新型馬氏體不銹鋼，利用Mo合金化提高N的固溶量并減少Cr、N化合引起的Cr消耗，進而提高耐蝕性[16]. 近些年隨著計算機技術(shù)的發(fā)展進步，高通量的材料計算和材料基因組的開發(fā)思路開始廣泛應(yīng)用于各種材料的研究中，例如具有合成或剝離可能性的2D材料[17]、高析氫反應(yīng)(HER)催化性能的二元合金材料[18]、高隔熱性能的ABO4白鎢礦[19]、高循環(huán)壽命的鋰離子電池電極材料[20]等等. 對于PEMFCs環(huán)境下的不銹鋼腐蝕性質(zhì)與設(shè)計改性相關(guān)的模擬計算研究，高通量方法會是高效可行的方案.

不銹鋼表相在PEMFCs的酸性環(huán)境下，表相鈍化層的成分以Cr2O3為主[21，22]，因此本文的計算以Cr2O3為建模基礎(chǔ)，對其進行了48種元素的摻雜替換，考慮了228種摻雜構(gòu)型、102種吸附情況、453種表面空位情況，并計算了摻雜后體系的的功函數(shù)、F原子吸附能、Cr/O空位形成能等腐蝕性質(zhì)，分析能夠提高不銹鋼鈍化層在PEMFCs環(huán)境耐蝕性的元素，以期為高耐蝕不銹鋼雙極板基材設(shè)計提供指導.

2 計算與建模方法

2.1 計算方法

本文采用基于密度泛函理論的第一性原理計算軟件VASP(Viennaabinitiosimulation package)[23]進行計算，贗勢采用平面投影綴加波(Project augmented wave，PAW)[24]，交換關(guān)聯(lián)能采用廣義梯度近似(Generalized gradient approximation，GGA)下的PBE(Perdew-Burke-Ernzerh)泛函[25]，K點采樣使用(4×2×1)的Monkhorst-Pack方法[26]，能量收斂設(shè)置為1×10-5eV，力收斂設(shè)置為1×10-2eV，截斷能設(shè)置為520 eV.

(1)替換能計算：本文中的摻雜形式主要考慮替換摻雜. 替換能是將材料中粒子替換為其他粒子的過程所需的能量，在本文中用于表征元素摻雜進入已有體系的難易程度，替換能越小表明摻雜越容易發(fā)生. 替換能的計算公式如下：

E(dis)=E(slab-X+Y)-

E(slab)-E(Y)+E(X)

(1)

式中E(dis)表示替換能，E(slab)表示未摻雜的切面模型能量，E(slab-X+Y)表示用Y元素替換切面模型中X元素后切面模型的能量，E(X)、E(Y)分別表示X、Y元素在各自穩(wěn)定體系中的單原子能量.

(2)功函數(shù)計算：功函數(shù)是電子從材料內(nèi)部移動到無限遠處所需能量，可以表明材料內(nèi)電子的穩(wěn)定程度，功函數(shù)越高體系越難失去電子. 功函數(shù)的計算公式如下[27]：

φ=Φ-Efermi

(2)

式中φ表示功函數(shù)，Φ表示真空電子能量，Efermi表示材料費米能級.

(3)吸附能計算：吸附能是吸附質(zhì)吸附到材料表面過程中需要的能量. 材料腐蝕模擬中計算腐蝕性粒子的在目標體系上的吸附能是研究目標體系在特定環(huán)境下耐蝕性的常用方法[28，29]，吸附能的大小可以表征體系受腐蝕性粒子影響的難易程度. PEMFCs環(huán)境下破壞鈍化層的腐蝕性粒子主要是F-[30，31]，因此本文采用F的吸附能表征摻雜體系耐受F-腐蝕的能力. 吸附能的計算公式如下[29]：

(3)

式中E(ads)表示吸附能，E(X+slab)表示吸附X原子之后切面模型的能量，E(slab)表示吸附之前的切面模型能量，E(X2)表示孤立X2分子的能量.

(4)空位形成能計算：不銹鋼PEMFCs雙極板的常見腐蝕失效形式是鈍化層溶解導致的點蝕[30，32]，也有文獻表明Cr2O3在酸性環(huán)境下且電極電位較高(約1.3V vs. SHE)時會出現(xiàn)溶解[21，33]. 材料溶解的難易程度常采用空位形成能來表征[34，35]，空位形成能是將材料表面移除一個原子形成空位所需要的能量. 空位形成能的計算公式如下[29]：

E(vac)=E(slab-X)+E(X)-E(slab)

(4)

式中E(vac)表示空位形成能，E(slab-X)表示移除X原子后切面模型的能量，E(X)表示穩(wěn)定的周期性體系中單個X原子的能量，E(slab)表示未移除原子時切面模型的能量.

2.2 建模方法

圖1 (a)Cr2O3的晶胞結(jié)構(gòu)；(b)Cr2O3的切面結(jié)構(gòu)(左視圖)Fig. 1 (a)Lattice structure of Cr2O3；(b)slab of Cr2O3 (left view)

圖2 (a)Cr2O3切面對稱性(俯視圖)；(b)選定周期的原子編號；(c)選定周期的吸附位點網(wǎng)格Fig. 2 (a)Symmetry of Cr2O3 slab (top view)；(b)atoms’ serial number of selected period；(c)Adsorption site grids of selected period

通過對研究周期中的4個Cr原子(Cr11，Cr12，Cr21，Cr22)進行替換實現(xiàn)摻雜，并通過計算替換能來判斷最容易發(fā)生摻雜的位點. 為了考察摻雜對于腐蝕性質(zhì)的影響，選擇最容易發(fā)生摻雜的位點作為摻雜位點，并且后續(xù)計算主要研究摻雜位點附近. 對于非金屬元素，額外計算了對周期內(nèi)O原子(O13，O14，O23，O24)替換摻雜的情況. 可以發(fā)現(xiàn)，Cr11、Cr22存在結(jié)構(gòu)上的對稱性，后續(xù)的計算也會驗證這一點；同樣(Cr12，Cr21)、(O11，O12)、(O21，O22)也是三組結(jié)構(gòu)對稱的位點.

吸附的初始位置設(shè)置在距離Cr2O3表面約2?處，并將研究周期劃分為如圖2c所示網(wǎng)格，計算每個格點的吸附位置和吸附能，篩選出多個未摻雜時容易發(fā)生吸附的位點，對于不同的摻雜位點選取最近的易吸附位點，計算摻雜對于吸附性質(zhì)的影響. 空位位點在摻雜位點周圍的Cr和O中選取，分別計算其空位形成能，以其中最低值作為摻雜點附近的空位形成能. 共計有48種元素參與了摻雜計算，將有害元素或者過于活潑的元素、惰性氣體元素等明顯難以實現(xiàn)摻雜的元素排除在外.

3 計算結(jié)果與討論

3.1 摻雜位點與類型

各元素摻雜的替換能如圖3所示，橫軸元素順序為元素周期表序列順序，替換能的大小體現(xiàn)摻雜發(fā)生的難易程度，替換能越低越容易發(fā)生摻雜，后續(xù)的各種耐蝕性參數(shù)的計算都基于最容易發(fā)生摻雜的位點進行. 可以發(fā)現(xiàn)，Cr11和Cr22位點的替換能幾乎完全相同，可以驗證這兩個位點在結(jié)構(gòu)上的對稱性，同樣Cr12和Cr21位點的替換能也幾乎相同；對于大部分金屬元素而言，都是Cr11/Cr22位點更容易發(fā)生摻雜，可能是因為Cr11/Cr22相對更接近材料表層. Mn(Cr11/Cr22：0.691 eV，Cr12/Cr21：0.649 eV)、Re(Cr11/Cr22：2.299 eV，Cr12/Cr21：2.092 eV)、Be(Cr11/Cr22：-1.339 eV，Cr12/Cr21：-1.414 eV)三種元素更容易發(fā)生Cr12/Cr21位點替換，因此這三種元素后續(xù)計算采用Cr12/Cr21位點摻雜. 另外可以看到Os、W、Tc的兩種Cr摻雜位點的替換能幾乎相同(兩種替換能差別<1%)，后續(xù)對于這些元素的計算將考慮兩種摻雜位點都可能發(fā)生的情況并取平均值.

圖3 Cr摻雜位點的替換能Fig. 3 Displacement energies of Cr doping sites

從過渡元素的Cr11/Cr22替換能可以發(fā)現(xiàn)較為明顯的周期性趨勢，大致以VI B族和VII B族為分界，前半部分元素Cr替換能為負，替換摻雜較易發(fā)生(Gd元素除外)；后半部分元素Cr替換能為正，替換摻雜較難發(fā)生.

對非金屬元素額外計算了替換O元素摻雜的替換能如圖4所示，C、N、P、S、Se、Te、F元素的O替換能比Cr替換能低約2～5 eV更容易發(fā)生對O的替換，針對這部分元素的后續(xù)計算采用O位點替換摻雜；B、Si兩種元素是金屬與非金屬之間的過渡性元素，相較于O替換更容易發(fā)生Cr的替換，后續(xù)計算采用Cr摻雜位點.

圖4 非金屬元素的Cr、O摻雜位點的替換能Fig. 4 Non-metal elements’ displacement energies of Cr/O

綜上所述，所有參與計算的元素可以分為Cr11/Cr22、Cr12/Cr21、O21/O22、O11/O12四類摻雜：Os、W、Tc為Cr11/Cr22或Cr12/Cr21摻雜；Mn、Re、Be為Cr12/Cr21摻雜；C、N為O21/O22摻雜；P、S、Se、Te、F為O11/O12摻雜；其余元素均為Cr11/Cr22摻雜.

3.2 耐蝕性參數(shù)計算結(jié)果

(1)摻雜后功函數(shù)計算結(jié)果：

各摻雜體系的功函數(shù)計算結(jié)果如圖5，未摻雜時的Cr2O3的功函數(shù)標注在圖中，其中22種元素摻雜后的功函數(shù)高于摻雜前. 由于Cr2O3本身的耐蝕性已經(jīng)比較優(yōu)秀，可以發(fā)現(xiàn)無摻雜的Cr2O3的功函數(shù)相對較高，摻雜后的元素中P、S、Co、Ni、Ga、Se、Ru、Rh、Te、Ir、Pt、Au的摻雜后功函數(shù)有較明顯的提升，C、N、Mg、V、Mn、Zn、Mo、Tc、Pd、In、W、OS的摻雜后元素與摻雜前變化不大，Cr2O3表相摻雜這些元素之后具有較高的穩(wěn)定性.

圖5 摻雜后的Cr2O3表面功函數(shù)Fig. 5 Work functions of Cr2O3 surface after doping (work function of Cr2O3 before doping is 4.7512 eV，as marked in the figure)

(2)F吸附能與摻雜后F吸附能變化量計算結(jié)果

在對研究周期劃分吸附位點網(wǎng)格之后，計算了網(wǎng)格中25個位點的吸附能，其中16個位點吸附于A吸附位點(圖6a中紫色位點)，吸附能為-3.572 eV；6個位點吸附于B吸附位點(圖6a中棕色位點)，吸附能為-3.573 eV；另外3個位點吸附于其他位置，吸附能為-1.804～-2.843 eV，因此這幾個位點較難發(fā)生吸附，后續(xù)只考慮A、B兩個吸附點(吸附至選定研究周期外的位點平移整數(shù)個周期至選定周期內(nèi)). A、B兩個位點的具體吸附位置如圖6b-e所示.

為研究摻雜對于吸附的影響，對于A、B兩點分別選取選擇與吸附位點距離最近的易摻雜位點進行替換摻雜并計算吸附能. 對于摻雜位點為Cr11/Cr22位點的元素，計算Cr11摻雜位點的A吸附位點吸附能(A點左移一個周期)和Cr22摻雜位點的B吸附位點吸附能；對于摻雜位點為Cr12/Cr21的元素，計算Cr21摻雜位點的A吸附位點吸附能和Cr12摻雜位點的B吸附位點吸附能；對于摻雜位點為O11/O12的元素，計算O12摻雜位點的A、B吸附位點吸附能(A點左移并上移一個周期)；對于摻雜位點為O21/O22的元素，計算O22摻雜位點的A吸附位點吸附能(A點上移一個周期)和O21摻雜位點的B吸附位點吸附能.

摻雜后的吸附能相較于未摻雜時的變化量如圖7，其中吸附能變化量為正表明相對摻雜前更難發(fā)生F的吸附，變化量為負表明F更容易吸附. 可以發(fā)現(xiàn)大部分元素兩個吸附位點的吸附能相差不大，只有摻雜類型為O11/O12過渡金屬元素摻雜的F吸附能變化量呈現(xiàn)出與Cr摻雜替換能、功函數(shù)相類似的周期性趨勢，大致以VII B族(Mn、Tc、Re)為界，前半部分F吸附能降低，后半部分反之. 金屬元素中，Co、Ni、Cu、Ge、Rh、Ag、In、Ir、Pt、Au對于F-的抑制作用較好；非金屬元素中的C、N、P、S、Se、Te均對F的吸附能有所提升，但是提升幅度不大.

圖 6 (a)選定研究周期內(nèi)的吸附位點分布；(b)A吸附點的吸附位置(俯視圖)；(c)A吸附點的吸附位置(側(cè)視圖)；(d)B吸附點的吸附位置(俯視圖)；(e)B吸附點的吸附位置(側(cè)視圖)Fig. 6 (a)Adsorption sites distribution of selected period (F atoms on purple/brown points adsorb to site A/B)；(b)adsorption position of site A (top view)；(c)adsorption position of site A (left view)；(d)adsorption position of site B (top view)；(e)adsorption position of site B (left view)

圖7 摻雜后的Cr2O3表面F吸附能變化量Fig. 7 F adsorption energy change of Cr2O3 surface after doping (compared to that before doping)

(3)Cr/O空位形成能與摻雜后Cr/O空位形成能變化量計算結(jié)果：

本文分別計算了摻雜對于Cr和O的空位形成能影響，未摻雜情況下研究周期內(nèi)的各位點空位形成能計算如表1，選擇最低值(即最容易形成空位)作為未摻雜時的空位形成能，即Cr空位形成能為0.917 eV，O空位形成能為4.126 eV.

表1 未摻雜的Cr2O3空位形成能

每種摻雜位點的元素均計算其周邊的Cr和O的空位形成能并取其中最低值作為空位形成能計算結(jié)果，選取范圍為距離摻雜位點3.5?以內(nèi)的表層原子. 不同摻雜類型元素采用的摻雜位點和空位位點篩選范圍如表2(Os、W、Tc分別考慮Cr11/Cr22和Cr12/Cr21兩種摻雜類型并取其均值，不作為單獨摻雜類型列出).

表2 各種摻雜類型選用的摻雜位點及其對應(yīng)的Cr、O空位篩選范圍

圖8 摻雜后Cr/O空位形成能變化量及各摻雜位點對應(yīng)的Cr/O空位位點Fig.8 Cr/O vacancy formation energy change after doping and corresponding doping sites of Cr/O vacancy (sites with lowest vacancy formation energy)

3.3 篩選與分析

摻雜后的F吸附能變化量、Cr空位形成能變化量、功函數(shù)的變化量的關(guān)系如圖9，以抑制F吸附和Cr空位形成、高穩(wěn)定性為目標，應(yīng)優(yōu)先考慮圖中右上區(qū)域中功函數(shù)較高的元素：C、N、In、Au、Ni、Rh、Pt、Co、Ir、Ru，計算結(jié)果表明這些元素在性能上可以滿足需求. 這些元素中，C、N、Ni是目前比較常見的不銹鋼雙極板的表面改性元素，也有不少學者對不銹鋼雙極板材料進行了離子注入[11，37]、滲C/N[38，39]等改性方法的探究，表明這三種元素可以提升不銹鋼表面的耐蝕性. Au、Rh、Pt、Ir則是典型的貴金屬元素，如果加入這些元素會導致制造成本過高，因此這四種元素的摻雜優(yōu)先級較低. In、Co、Ru三種元素具有較好耐蝕性參數(shù)，且價格適中，具備摻雜的可行性，這三種元素目前較少有應(yīng)用到不銹鋼雙極板的改性中；關(guān)于這三種元素對于不銹鋼材料的作用，In2O3常被用來作為不銹鋼的陰極光電化學保護材料[40]，Co的合金涂層可以提高不銹鋼材料的耐熱酸蝕性能[41]，Ru的加入則可以使不銹鋼的耐蝕性能和抗應(yīng)力腐蝕開裂性能得到提高[42]，因此In、Co、Ru是值得嘗試的不銹鋼雙極板改性元素；另外Co、Ru的最低Cr摻雜替換能偏高，實現(xiàn)摻雜的條件可能會更苛刻. 這些元素中，有文獻指出Ni2+、Co2+會毒化質(zhì)子交換膜中的催化劑，導致電池性能下降[5]，因此，如果Ni、Co的摻雜難以將表面的耐蝕性提升至幾乎不會有離子析出的程度，則應(yīng)謹慎考慮選用此兩種元素.

圖9 摻雜后F吸附能變化量、Cr空位形成能變化量和功函數(shù)的關(guān)系Fig.9 The relationships among F adsorption energy change，Cr vacancy formation energy change and work function after doping

除了上述綜合性能都較好的元素，也存在不少部分耐蝕性參數(shù)較好而某項參數(shù)不高的元素，例如Ag、Cu、Sn、Ge，除了摻雜后的功函數(shù)稍有下降外，Cr空位形成和F吸附的抑制效果均較好，其中對于Ag有學者通過離子注入Ag的方式對不銹鋼雙極板表面改性，使耐蝕性得到提升[12]. Mg、P、S、Zn摻雜對于F吸附能有一定提升且功函數(shù)較高，而對于Cr空位形成能的改善效果一般，其中P、S是不銹鋼中典型的有害元素，且可能會對質(zhì)子交換膜中的催化劑造成不利影響[5]，選用優(yōu)先級較低；Mg2+同樣會對質(zhì)子交換膜中的催化劑造成破壞，選用優(yōu)先級較低.

為探究不同元素摻雜后耐蝕性變化規(guī)律的原因，對部分元素摻雜后的材料表面電子態(tài)密度(Density of states，DOS)進行了計算. 對于金屬元素摻雜，耐蝕性改善效果較好的Ni和較差的Zr元素電子態(tài)密度如圖10a-b，摻雜的Ni的電子主要分布于費米能級以下約2 eV至費米能級處，Ni的d電子軌道與Cr2O3中Cr、O原子的電子軌道產(chǎn)生雜化，表明Ni與Cr2O3基體中的原子形成了較強的鍵合作用[27]，進而使得Cr脫離基體形成空位所需能量增加，也更難與F形成新的鍵合也即F的化學吸附更難以發(fā)生；對于摻雜后耐蝕性參數(shù)降低的Zr，可以發(fā)現(xiàn)相比于Ni的電子能態(tài)分布，Zr的d軌道電子主要分布于費米能級以上，s、p軌道電子則分別在約-50 eV和-28 eV處呈現(xiàn)出明顯的局域化現(xiàn)象，表明Zr與Cr2O3中的原子未能形成較強的鍵合作用，進而導致了Cr空位、F吸附更容易發(fā)生，同時d軌道電子在費米能級以上的大量分布也表明相較于Ni摻雜，Zr摻雜的電子更容易逸出，這也解釋了Zr摻雜后功函數(shù)的降低.

圖10 摻雜后Cr2O3表面的電子態(tài)密度：(a)Ni摻雜；(b)Zr摻雜；(c)C摻雜；(d)Si摻雜Fig. 10 DOSs of Cr2O3 surface after being doped：(a)DOS of Ni-doped surface；(b)DOS of Zr-doped surface；(c)DOS of C-doped surface；(d)DOS of Si-doped surface

非金屬元素中耐蝕性改善效果較好的摻雜元素C和較差的Si的電子態(tài)密度如圖10c-d，與金屬元素摻雜不同的是，C、Si的電子態(tài)密度分布表明這兩種摻雜元素的原子軌道均與基體原子軌道發(fā)生雜化，其中C的態(tài)密度峰趨于與Cr的態(tài)密度峰重合，而Si的態(tài)密度峰趨于與O的態(tài)密度峰重合，s軌道的分布尤為明顯. 這一現(xiàn)象主要由不同的摻雜類型導致，C采用的摻雜類型替換的是O原子，則C會趨于與原本與O成鍵的Cr原子成鍵；Si則是采用了Cr原子替換摻雜，則Si會趨于與被替換的Cr周圍的O原子成鍵，這一成鍵趨勢的差異導致了C和Si的引入對于Cr空位、F吸附的不同影響.

4 結(jié) 論

(2)綜合比較各種腐蝕參數(shù)，C、N、In、Ru摻雜可以提高功函數(shù)、抑制F吸附和Cr空位形成，且對于質(zhì)子交換膜的潛在危害較低；Au、Rh、Ir、Pt、Ni、Co摻雜的耐蝕性參數(shù)計算結(jié)果較好，但是存在價格、危害質(zhì)子膜等問題. Ag、Cu、Sn、Ge摻雜對Cr空位形成和F吸附的抑制效果均較好，功函數(shù)稍有下降；Zn摻雜后的功函數(shù)數(shù)值較好且有利于抑制F吸附.

(3)從電子態(tài)密度計算結(jié)果推斷摻雜原子與基體原子之間的成鍵作用會對摻雜后的耐蝕性產(chǎn)生影響，摻雜原子與Cr之間形成穩(wěn)定鍵合是摻雜原子提升耐蝕性的重要原因.