蔣紅星 孫俁哲 孫 升 張統(tǒng)一

（上海大學材料基因組工程研究院，上海 200444）

20世紀以來，金屬材料的使用量急劇增長，其開采冶煉加工過程消耗了大量的能源，每年因為電化學腐蝕造成了大量金屬材料的浪費，所以金屬的腐蝕溶解問題得到了廣泛的關注。金屬電極的腐蝕受應力應變和電化學環(huán)境的影響，因此有必要研究應變-電荷耦合場下金屬電極的腐蝕特性。Gutman等［1］研究發(fā)現(xiàn)，外加力場會導致金屬材料電化學勢的改變，在金屬電極發(fā)生彈性變形時，隨著拉伸／壓縮應力的增加，金屬的開路電壓呈線性降低的趨勢。這是由于力學化學作用使應變能轉變?yōu)殡娀瘜W能，即金屬的表面活性增大，腐蝕傾向也增大。國內關于應變腐蝕的研究報道［2-3］也很多，如 Zhang等［3］研究了拉應力作用下X80鋼在酸性溶液中的電化學行為，發(fā)現(xiàn)隨著拉應力的增加，電極的開路電位降低，反應活性增加。Mavrikakis等［4］通過密度泛函理論（density functional theory，DFT）研究了應變對金屬電極表面反應活性的影響，計算中對Ru（0001）晶面施加面內應變，發(fā)現(xiàn)表面應變與O原子／CO分子的吸附能和活化能壘之間存在關聯(lián)性，相對于無應變表面，拉伸應變表面降低了對O原子／CO分子的吸附能以及CO分解過程的能壘，并將這種效應歸因于金屬d能帶中心隨應變的變化。

鋅具有蘊藏量豐富、低毒、無污染、高能量密度以及低成本等特點，而被廣泛用作電池負極材料。其中可充電鋅離子電池（zinc-ion battery，ZIB）由于其良好的環(huán)境兼容性、極高容量密度、低廉的價格和很高的安全性，被認為是今后大規(guī)模能量儲存的方向之一［5-6］。鍍鋅是一種犧牲陽極保護陰極的電化學防腐蝕方法，其工藝成熟、成本較低，被廣泛應用于制造業(yè)［7］。

但目前的試驗方法無法從微觀層面研究腐蝕過程中原子的溶解和電荷的轉移過程，而跨尺度計算則可以通過原子計算結合宏觀溶液計算較好地模擬鋅的電化學溶解。Schwarz等［8］通過使用耦合密度泛函理論／連續(xù)介質溶液方法，在寬電壓范圍內的密集計算，揭示了吸附在金屬晶體臺階邊緣上水的部分氧化機制。本文選擇密度泛函／連續(xù)介質溶液方法，對鋅電極的溶解腐蝕過程進行了詳細研究，并結合符號回歸算法預測了鋅電極發(fā)生溶解腐蝕的臨界電極表面應變-電荷組合。

1 數(shù)學建模與方法

圖1 電極／溶液模型和鋅電極／溶質鋅原子結構示意圖Fig.1 Schematic of the electrode／solution model and zinc electrode／solute zinc atoms

計算系統(tǒng)在整個電化學回路中的位置如圖1所示，電極為Zn，并浸入含NaF的電解質溶液中。電極使用DFT計算，電解質溶液采用連續(xù)介質模型計算［9］，計算采用超胞方法。電極模型采用了22個Zn原子，其中20個原子形成一個含有5層原子的平板電極（每層有2×2個Zn原子），另2個鋅溶質原子分布在平板的兩側，形成對稱結構，如圖1中的插圖所示。Zn電極與電解質溶液接觸的表面方向為密排六方晶格的＜0001＞方向，［0001］方向被定義為 z軸，Zn電極最外層的原子面則位于z＝0。電極兩側填充電解質溶液，超胞沿z向的長度被設定為44.8?，電極厚度為10.28?。

本文使用開源軟件 JDFTx［10-13］進行密度泛函／連續(xù)介質溶液耦合計算，電解質溶液模型采用SaLSA模型，該模型具有非局域介電響應特性［9］，對電解質溶液中強極化和帶電溶質的計算的準確度較高。使用濃度為1 M的隱性Na+和F-離子為鋅電極提供徳拜屏蔽效應。DFT計算采用平面波展開，平面波矢空間的K點網(wǎng)格劃分為12×12×2，截斷能量設定為707 eV，并采用周期性邊界條件，交換關聯(lián)函數(shù)使用廣義梯度近似（general gradient approximation，GGA）［14］，鋅 原 子 使 用Ultrasoft贗勢方法［15］。電極的電子態(tài)密度計算使用線性四面體插值方法［16］。

計算前先進行結構優(yōu)化，計算步驟如下：（1）將建好的Zn電極模型去除2個溶質原子后放入溶液中進行晶格參數(shù)優(yōu)化（即整個結構處于能量最小化），得到Zn電極HCP結構的平均晶格參數(shù)為：；（2）將兩個溶質Zn原子放入第1步的結構兩側進行能量最小化（第1步中的結構保持固定，僅松弛2個溶質原子的位置），溶質原子會吸附在鋅電極表面3個相鄰Zn原子的幾何中心上方。在之后一系列計算中電極的晶格參數(shù)和溶質原子在平行于電極表面內的位置都保持不變，改變溶質Zn原子與電極表面的距離（從 z＝1．55?至 z＝14．65?）進行一系列能量計算。在密度泛函計算中，金屬表面的應變通常是通過改變表面方向的晶格參數(shù)來實現(xiàn)的［4，17-19］，本文也采用該方法，即通過同時改變晶格a和b方向的參數(shù)來施加面內應變。

DFT計算時將凈電荷施加在鋅電極上，因此凈電荷位于電極和2個孤立Zn原子中。表面電荷由電極最外表面原子的電荷密度相對于電極內部原子的電荷密度的變化來定義，通常使用吉布斯突變表面模型，定義為表面每單位表面積的凈電荷。為了區(qū)分這兩類凈電荷，將整個計算系統(tǒng)的凈電荷簡稱為凈電荷，而表面的凈電荷稱為表面電荷。凈電荷是由施加的電位引起的，而表面電荷是由自由表面引起的。表面電荷隨凈電荷的變化而變化。

密度泛函／連續(xù)介質溶液計算中分別使用自洽場（SCF）和Polak-Ribière共軛梯度法進行電子和電解質的弛豫。在兩次連續(xù)的SCF自洽循環(huán)中，能量收斂的相對準確度被設定為10-7Hartree（1 Hartree≈27.211 4 eV）。采用 L?wdin電荷布居分析［21-22］計算電極和溶質Zn原子的電荷分布。

由于計算Zn溶解／吸附的臨界（凈電荷、電極面內應變）組合需花費較多時間和成本，因此基于機器學習中的符號回歸算法（symbolic regression）［23-24］，找出擬合訓練數(shù)據(jù)的函數(shù)關系式，然后利用函數(shù)關系式，求得具有給定性質的參數(shù)組合，進而加速尋找臨界凈電荷-電極表面應變組合。

2 結果和討論

為了計算方便，采用22個鋅原子的凈電荷總值表征系統(tǒng)的凈電荷。圖2（a）和圖2（b）分別是在凈電荷為3.5｜e｜和4.0｜e｜時，不同電極面內應變下能量隨電極-溶質原子距離的變化（圖中實線表示使鋅溶質原子自動溶解的應變）。能量Ead的定義為：

式中：E［Zn（slab）+Zn（z→∞）］表示計算系統(tǒng)的能量，此時溶質鋅原子遠離電極表面，Zn（z→∞）］由 Zn（z＝14.7?）近似，在 z≥7?時溶質鋅原子和Zn電極之間的相互作用可忽略不計。

如圖2（a）所示，在凈電荷為 3.5｜e｜、未施加應變的條件下，在2.46?＜z＜4.5?范圍內，系統(tǒng)能量處于1 eV的平臺，即在此范圍內溶質原子的溶解／吸附近似處于平衡狀態(tài)。在電極面內施加5%拉伸應變時，在1.6?＜z＜4.5?范圍內出現(xiàn)了一段能量峰值，即從2.46?處的能量低點（電極表面的吸附位置）到4.5?處的能量高點，約0.17 eV的能量勢壘，該勢壘會阻礙溶質原子溶解進入電解質溶液。在電極面內施加 -5% 應變（壓縮應變）時，吸附能單調遞減，說明表面溶質原子產生了自發(fā)溶解。z＞4.5?時，不同應變下系統(tǒng)能量均不斷降低，當z＞6.9?時基本穩(wěn)定，表明溶質原子在遠離電極表面時與電極間的相互作用微弱，而在4.5?＜z＜6.9?時能量明顯降低，說明溶液中的溶質原子將難以跨越能壘而傳遞到電極上。如上所述，（3.5｜e｜，0）為一組臨界組合（電荷，應變）。

如圖2（b）所示，在凈電荷為 4｜e｜、未施加應變的條件下，系統(tǒng)能量隨溶質原子-電極表面間距增加而單調減小，說明表面溶質原子產生了自發(fā)溶解，且這一過程不可逆。在電極面內施加5%和15%拉伸應變時，在2.46?＜z＜4.5?范圍內，能量下降坡度隨著拉伸應變的增加逐漸減緩，此處能量坡度定義為溶質原子從z＝4.5?進入到吸附位置z＝2.46?時的能量勢壘，該值在0%、5%和15%應變下分別為1.02、0.63和0.11 eV。在施加20%拉伸應變時，能量坡度為負值，在z＝2.46?時最低為0.61 eV，在z＝4.5?時最高為0.82 eV，能量勢壘為0.21 eV，阻礙溶質原子進入電解質溶液。在z＞4.5?范圍內，能量變化趨勢與圖2（a）的相似，4.5?＜z＜6.9?范圍內的能量坡度阻止了溶質原子吸附到電極表面。

圖2 不同凈電荷下能量隨電極面內應變的變化Fig.2 Energy profiles of zinc atoms on a zinc electrode under different surface strains at different net charges

圖3（a、b）和圖 3（c、d）分別是在凈電荷為3.5｜e｜和 4.0｜e｜時，不同電極面內應變下溶質 Zn原子和電極最外表面鋅原子所帶電荷隨電極-溶質原子距離的變化。圖3（a）中，溶質原子的最大電荷均出現(xiàn)在z＝4.4?時，在0%、5%、-5%應變下該數(shù)值分別為 1.65｜e｜、1.53｜e｜和1.76｜e｜。在［1.6，4.4］?范圍內，不同應變下每個溶質Zn原子所帶電荷均隨著電極-溶質原子間距的增加而增加，如在5%拉伸應變下，溶質原子所帶電荷從1.11｜e｜增加到1.53｜e｜。圖3（b）中，在［1.6，4.4］?距離范圍內，不同應變下電極最外表面每個鋅原子所帶電荷均隨著電極-溶質原子間距的增加單調減少，在0%應變下，鋅原子所帶電荷從0.32｜e｜減小到 0.07｜e｜，5%應變下從0.19｜e｜減小到 0.08｜e｜，在 -5%應變下從 0.16｜e｜減小到0.03｜e｜。當z＞6.9?時，在壓縮應變和無應變條件下，電極最外表面鋅原子所帶電荷接近，約0.09｜e｜，在拉伸應變下，約 0.11｜e｜。

圖3（c）中，隨著拉伸應變的增加，溶質鋅原子所帶電荷降低，如在z＝6.9?，電極面內拉伸應變?yōu)?%、5%、15%和20%時，溶質鋅原子所帶電荷分別為 1.89｜e｜、1.81｜e｜、1.63｜e｜和 1.53｜e｜。圖3（d）中，在［1.6，3.4］?范圍內，0%應變下，電極最外表面鋅原子所帶電荷最高，如在z＝2．56 ?時為0.23｜e｜，5%、15%、20%應變下分別對應為0.17｜e｜、0.18｜e｜、0.19｜e｜，該結果與圖 3（b）中相似，即應變導致在［1.6，3.4］?范圍內電極最外表面原子所帶電荷下降，但拉伸應變的增加導致電荷升高，壓縮應變的增加導致電荷下降。在z＞4．4?時，隨著拉伸應變的增加，最外表面原子所帶電荷升高且隨著溶質鋅原子與電極表面距離的增加基本不變，如在 z＝6.9?時，0%、5%、15% 和 20%下分別為 0.07｜e｜、0.10｜e｜、0.12｜e｜、0.13｜e｜。結合圖2與圖3（a）和3（c）中的實線部分，可見在3.5｜e｜凈電荷和-5%應變、4.0｜e｜凈電荷和 0%應變、4.0｜e｜凈電荷和 5%應變、4.0｜e｜凈電荷和15%應變條件下，溶質鋅原子溶解進入電解質溶液所帶電荷分別為1.45｜e｜、1.89｜e｜、1.81｜e｜、1.63｜e｜，比較能量曲線發(fā)現(xiàn)，溶質原子所帶電荷越接近經(jīng)典理論的+2價態(tài)，溶質Zn原子越易溶解。

為了尋找凈電荷4｜e｜下溶質鋅原子溶解的臨界應變值，使用SR進行符號回歸訓練，圖4顯示了在凈電荷為4｜e｜時，不同溶質原子-電極表面距離和拉伸應變條件下系統(tǒng)能量的符號回歸數(shù)值模型。將訓練過程中的4組數(shù)據(jù)（共40個，圖4中黑色數(shù)據(jù)點）均作為訓練集，每個數(shù)據(jù)包含了溶質原子-電極表面間距與電極面內應變2個特征量和1個系統(tǒng)能量輸出量。使用 gplearn（https：／／gplearn.readthedocs.io）中的符號回歸模塊對歸一化后的數(shù)據(jù)集進行訓練，符號回歸模型使用的參數(shù)有：population_size＝5000，迭代次數(shù)（generations）為 113，計算精度（stopping_criteria）設置為0.008，交叉概率（p_crossover）為0.7，子樹變異概率p_subtree_mutation＝0.1，hoist變異概率p_hoist_mutation＝0.05，點變異概率p_point_mutation＝0.1，最大樣本抽取 max_samples＝0.9，過程打印 verbose＝1，稀疏系數(shù) parsimony_coefficient＝0.000 12。根據(jù)圖4中符號回歸訓練模型，在［0%，20%］應變范圍內等間距取100個應變值并代入圖4的符號回歸模型，比較發(fā)現(xiàn)在應變?yōu)?6.75%時，溶質原子在電極表面的溶解／吸附最接近于臨界狀態(tài)。

圖3 在凈電荷為3.5｜e｜和4.0｜e｜時，不同電極面內應變下溶質Zn原子和電極最外表面鋅原子所帶電荷隨電極-溶質原子距離的變化Fig.3 Charge curves of solute zinc atoms and zinc atoms on the outermost surface of electrode as a function of distance from electrode to solute atoms under different electrode surface strains when the net charges are 3.5｜e｜and 4.0｜e｜

圖4 凈電荷為4｜e｜時的符號回歸訓練模型Fig.4 Symbolic regression trainingmodel under net charge of 4｜e｜

為了驗證SR預測的臨界組合（電荷，應變）的準確性，進行了凈電荷為4｜e｜、電極面內應變?yōu)?6.75%條件下的計算。如圖5所示，實線為符號回歸模型預測值，方塊為采用JDFTx計算的電極表面應變?yōu)?6.75%時的數(shù)值，可見計算值與符號回歸模型預測值基本吻合，說明了符號回歸模型的預測精度較高。在2.46?＜z＜4.5?范圍內，模型預測值和計算值都為0.61 eV左右，沒有出現(xiàn)能量勢壘，即溶質原子在吸附位置處于臨界狀態(tài)。在z＞6.9?時，兩者均逐漸穩(wěn)定，說明符號回歸模型可以預測溶質原子進入電解質液后與電極表面微弱的相互作用。結合JDFTx計算和符號回歸的預測值，得到了兩組溶質原子溶解／吸附臨界組合（電荷，應變），分別為（3.5｜e｜，0%）和（4｜e｜，16.75%）。從圖中還能看出，當電荷小于6.9?時，另外兩組以電極表面應變?yōu)樽宰兞康那€的應變值存在差異，而在z＞6.9?后逐漸與電極表面應變16.75%時的曲線重合，并趨于穩(wěn)定，與上述解釋相符。

圖5 符號回歸算法預測曲線和JDFTx計算值Fig.5 Curves predicted by using symbol regression algorithm and values calculated by JDFTx

3 結束語

本文使用密度泛函／連續(xù)介質溶液模型耦合計算研究了應變-電荷耦合場下Zn電極的溶解機制。在凈電荷3.5｜e｜下分別在電極表面內施加了5%（拉伸）和-5%（壓縮）應變，鋅電極表面溶質原子的近似溶解／吸附平衡狀態(tài)被打破，施加5%應變抑制了溶質Zn原子的溶解，而施加-5%應變則促進了溶質Zn原子的溶解。此外，在系統(tǒng)凈電荷為4.0｜e｜以及施加不同拉伸應變的條件下，拉伸應變越大抑制作用越明顯，直至出現(xiàn)能量勢壘阻止了溶質Zn原子的溶解，且拉伸應變的增加導致溶質Zn原子進入電解質溶液所帶電荷也下降。結合JDFTx計算值和符號回歸模型，可以預測固定凈電荷下Zn電極表面溶質原子溶解／吸附平衡狀態(tài)的臨界應變。但是本文計算得出的拉伸應變抑制溶解的結果與通常認知不符，可能的原因是通常試驗采用的是宏觀多晶材料，而本文計算的是一種完美單晶金屬材料。