金屬腐蝕的多尺度計算模擬研究進展

魏 薪，董超芳，徐奧妮，李 妮，肖 葵，李曉剛

(北京科技大學 腐蝕與防護中心 教育部腐蝕與防護重點實驗室，北京 100083)

1 前 言

金屬的腐蝕受到材料因素、環境因素和界面狀態的影響，任何單一因素或耦合條件的改變都可能造成腐蝕規律的變化[1]，影響因素的復雜程度造成了腐蝕行為的多樣性，其機理研究的難度也正源于此。伴隨著計算機的發展，“計算和模擬”已經成為了繼實驗、形式理論之后的第3條發現新概念和機理的途徑[2, 3]。在通過實驗建立的腐蝕模型基礎上，可以實現對點蝕[4-6]、晶間腐蝕[7]、應力腐蝕[8]、沖蝕[9, 10]、擴散傳質[11, 12]和電極過程動力學[13]等的理論研究。近年來，隨著試驗觀測技術的提高，掃描電化學顯微鏡(SECM)[14]、原子力顯微鏡(AFM)[15]、掃描開爾文力顯微鏡(SKPFM)[16]等微觀觀測儀器已經應用到了腐蝕研究中。這些儀器的觀測尺度一般在1 mm到1 μm的范圍內，可以細致到觀測單個蝕坑的形貌，研究單個蝕坑周圍的電位，這對揭示腐蝕的深層機理有很大的促進。觀測手段的發展，為介觀尺度下的腐蝕過程的建模和仿真提供了試驗數據基礎。

另一方面，受反應時間和尺度限制，大多數腐蝕實驗研究方法，并不能直接觀察材料的腐蝕萌生行為。研究者只能通過實驗獲得的數據進行分析進而推斷反應發生過程，這一方法容易發生錯誤致使人們不能正確獲得真正的腐蝕機理，從而不能正確指導材料的設計與使用。原位測試技術的發展在某種程度上解決了部分問題，但由于其對試樣和實驗條件要求極高并不能廣泛的應用。近年來，理論計算用軟件得到了迅速發展，使得腐蝕機理的研究取得了突破性的進展，對材料的設計和使用具有積極的指導作用。目前，腐蝕領域中應用較廣泛的計算方法包括基于密度泛函的第一性原理和分子動力學方法、基于牛頓力學的分子動力學計算、基于電化學的熱力學和動力學計算以及基于斷裂力學的有限元方法等。

計算方法的選擇與所研究對象的尺度有關，已有從原子到宏觀尺度的方法包括：第一性原理計算、Monte Carlo方法、分子動力學、元胞自動機、有限元、邊界元等，其中第一性原理計算方法受到較大關注。Di Stefano等[17]利用基于密度泛函的第一性原理計算方法研究了∑3(111)[-110]晶界和∑5(210)[001]晶界與H的相互作用，結果表明：這兩種晶界與H原子的相互作用行為截然不同。密排的∑3晶界既不能捕獲H原子也不能促進其在晶界的擴散，在材料內部充當二維的擴散阻擋層；而∑5晶界能夠提供大量的H溶解位，且相比于塊體材料，H原子沿∑5晶界的擴散系數提高了兩個數量級。Kulkoua等[18]利用第一性原理計算方法研究了H與B2結構TiFe合金的相互作用，通過計算H原子晶界和自由表面的吸附能和偏聚能表明，H有助于晶間開裂的發生。

本文將對不同尺度下金屬腐蝕研究的計算方法進行總結，并以大氣環境下鋁合金的腐蝕與氫致開裂為計算對象，利用第一性原理和內聚有限元的跨尺度計算方法，研究了氫原子偏聚導致的晶間裂紋的萌生與擴展行為，將第一性原理計算獲得的晶界結合能輸入到內聚有限元計算中，實現了從原子尺度到宏觀尺度的跨尺度計算。

2 計算方法在金屬腐蝕研究中的應用

2.1 第一性原理計算方法及其在腐蝕研究中的應用

隨著量子理論的建立和計算機技術的發展，人們希望能夠借助計算機對微觀體系的量子力學方程進行數值求解，然而量子力學的基本方程—薛定諤方程的求解是極其復雜的。然而，電子密度泛函理論(DFT)的確立[19, 20]克服了這種復雜性。電子密度泛函理論是20世紀60年代在Thomas-Fermi理論的基礎上發展起來的量子理論的一種表述方。傳統的量子理論將波函數作為體系的基本物理量，而密度泛函理論則通過粒子密度來描述體系基態的物理性質。密度泛函理論是一種完全基于量子力學的從頭算理論，但是為了與其他的量子化學從頭算方法區分，人們通常把基于密度泛函理論的計算叫做第一性原理。

在密度泛函理論中，所有的近似都被集中到稱為交換相關能的一項上，所以密度泛函理論的精度直接由交換相關能量泛函的近似形式決定。尋找更好的交換相關近似就成為密度泛函理論體系發展的一條主線。交換相關能量泛函的一個最初的簡單近似方法是局域密度近似(LDA)，即用具有相同密度的均勻電子氣的交換相關泛函作為對應的非均勻系統的近似值，這樣一個簡單近似往往能給出很好的結果。如果進一步分別考慮不同自旋分量的電子密度，可得到自旋極化的局域密度近似L(S)DA。盡管L(S)DA獲得了巨大的成功，但是也有許多不足之處，比如系統地高估結合能。在L(S)DA基礎上的改進有廣義梯度近似(GGA)。在GGA近似下，交換相關能是電子密度及其梯度的泛函。構造GGA交換相關泛函的方法分為兩個流派。一個是以Becke 為首的一派，認為“一切都是合法的”，人們可以任何理由選擇任何可能的泛函形式，而這種形式的好壞由實際計算來決定。通常，這類泛函的參數由擬合大量的計算數據得到。另外一個流派以Perdew 為首，他們認為發展交換相關泛函必須以一定的物理規律為基礎，這些規律包括標度關系、漸進行為等。基于這種理念構造的一個著名的GGA泛函是PBE泛函[21]，也是現在用得最廣泛的GGA泛函之一。電荷密度泛函理論在局域電荷密度近似或是梯度修正的版本，這是由Perdew和Wang所發展的GGA。DFT所描述的電子氣體交互作用被認為是對大部分的狀況都是精確的，并且他是唯一能實際有效分析周期性系統的理論方法。梯度修正的方法在研究表面的過程，小分子的性質，氫鍵晶體以及內部空間的晶體是比較準確的。眾所皆知，LDA會低估分子的鍵長或者鍵能，以及晶體的晶格參數，而GGA會補救這些缺點。而GGA在里斯晶體會過度修正LDA結果。當LDA與實驗符合非常好的時候，GGA會高估晶格長度。

利用第一性原理計算方法研究金屬材料在原子尺度的腐蝕行為對探究腐蝕機理有著十分重要的意義。基于密度泛函的計算方法可從原子之間相互作用的角度計算金屬材料的固有性質，或描述金屬基體或鈍化膜表面與環境中分子或離子等之間的化學反應。

2.1.1 表面吸附導致鈍化膜的形成與破壞

金屬暴露于各類環境下，與氣體或溶液中的氣體和離子發生相互作用，具體宏觀表現為金屬氧化膜的形成與破壞以及金屬表面與腐蝕性氣體和離子的相互作用。由此可見，吸附對腐蝕行為影響顯著，金屬表面對氣體或離子的吸附分為物理吸附和化學吸附。判斷被吸附物與表面是否發生化學反應需要計算體系的電子性質，如態密度、電荷密度、電子局域函數、布局分析、Bader電荷分析等。

本課題組[22]利用密度泛函方法計算了H2O和SO2在面心立方金屬Cu(111) 表面的單吸附和共吸附行為。通過對幾何構型，能量，電荷轉移及態密度的分析可得到以下結論。單分子H2O和SO2均不能與化學鍵的形式吸附在表面。對于共吸附構型，當覆蓋度為0.25 ML時，H2O和SO2分子分別吸附在Cu(111)表面；當覆蓋度增大到0.5 ML時，H2O分子解離為OH和H，其中OH吸附在表面，而H與SO2鍵合遠離表面。

Guo等[23]利用第一性原理計算表明，H2O分子在Al(111)純表面的解離能為248.32 kJ/mol，然而當存在預吸附O時，解離能降低到128.53 kJ/mol，說明O原子的存在大大降低了H2O在鋁表面的解離能。基于密度泛函的計算能夠有效的研究H2O與金屬的相互作用，但是在計算中卻存在一個嚴重的問題，那就是標準的密度泛函方法無法描述H2O分子與金屬表面的弱相互作用[24]。目前，已有大量研究[25,26]證明在考慮H2O分子與金屬表面相互作用時應考慮范德華色散力的作用。本課題組[27,28]采用范德華密度泛函計算方法研究了O2和H2O分子在Al(111)表面的單吸附和共吸附行為。結果表明，O2分子可自發解離，以O原子的形式吸附在表面，伴隨著較大的吸附能；而H2O分子只能以分子的形式弱吸附在表面。同樣地，預吸附的O原子僅能促進H2O分子的變形，大大降低其解離能壘，不能使其解離吸附。當O2分子存在時，且兩個O原子同時與H2O分子中的H原子產生相互作用時，H2O分子解離為OH和H，吸附在金屬Al(111)表面。通過對O2和H2O分子的共吸附行為的研究可知，Al表面鈍化膜的形成機理為H2O分子的存在促進了吸附的原子O由表面向內部的滲透而增厚，而不是金屬原子向外部的遷移。

Digne等[29]利用第一性原理計算的方法建立了γ-Al2O3表面模型。利用密度泛函方法，作者詳細研究了γ-Al2O3三種表面(100)、(110)和(111)的酸基性質，并考慮了與溫度相關的羥基化表面的覆蓋度。模擬與大量實驗結合解決了在紅外光譜中出現的OH-伸縮頻率的問題。研究還確定了酸化表面位的精確本性。酸化的強度可以通過模擬目標分子，如CO，和計算表面電子性質來量化。這些結果嚴重挑戰了以前的無缺陷尖晶石γ-Al2O3的表面模型并建立了γ-Al2O3酸化性質更為嚴格的描述。本課題組的研究成果[30]利用考慮強關聯效應的密度泛函方法研究了SO2分子在金屬鎳表面鈍化膜NiO(111)表面的吸附行為，結果表明，表面O空位的存在大大加強了表面對SO2分子的吸附，進一步說明表面鈍化膜的完整性對于材料耐蝕性的關鍵作用。這一系列的研究充分說明第一性原理可用于鋁合金與大氣環境中氣體相互作用的計算。

2.1.2 氫原子在界面處的偏聚

材料內部存在大量的晶界以及第二相與基體的界面，這些界面屬于高能區。在腐蝕環境下，腐蝕優先發生在界面處，且氫原子在界面的聚集導致氫脆的發生。因此，計算界面的結合能以及氫原子聚集對界面結合能的影響對研究應力腐蝕、氫致開裂等局部腐蝕機理具有重要的意義。Fernandea等[31]利用密度泛函與基于過渡態理論和熱力學的靜態模型相結合的方法研究了氫原子與金屬鎢的相互作用。該模型得到的與溫度相關的數據有助于理解宏觀實驗獲得的結果。研究結果有效的解釋了實驗測得的H擴散系數與密度泛函計算獲得的結果存在偏差是由于空位濃度與溫度的差異。Kristinsdóttir等[32]利用密度泛函方法研究了氫原子在23種過渡態金屬密排表面的擴散行為。作者選擇了第4、5和6周期，從Sc到Au，除了Mn、Tc和Hf的d-金屬。氫原子在這些金屬表面的勢能面被建立，且從一個極小值到另一個極小值的能壘與微動彈性帶計算相比較。結果表明，所有的擴散激活能相對低或從Pt/0.04 eV到Y和Zr/0.28 eV。作者還估計了隧道效應開始發生的溫度。Ferrina等[33]利用周期自洽DFT-GGA(PW91)方法研究氫原子與17個過渡族金屬表面，即面心立方金屬的(100)和(111)面，密排六方金屬的(0001)和體心立方金屬的(100)和(110)的相互作用。計算得到與實驗結果一致的幾何結構和結合能。對于同一金屬的密排面和較開放的表面，計算得到的結合能有很大的不同。同一金屬不同亞表面的氫原子幾何結構基本一樣，但是結合能有很大不同。除了Ag和Au表面H的吸附相對于H2時放熱的。對于所有研究的金屬，H在亞表面層的優先吸附位沒有在表面的優先吸附位穩定。H從表面到第一層亞表面的擴散激活能由兩個表面熱力學穩定性的差異控制。擴散從第一亞表面擴散進入下一層通常沒有很大的熱力學能壘，但有適當的動力學能壘。對于同一種金屬，H在不同表面向內擴散的激活能壘也有很大的不同。Mao等[34]利用第一性原理和分子動力學方法研究了間隙原子H在Er2O3內的擴散，計算了H在塊體材料，(001)表面和 Σ13 (4-3-1)/[111] 扭轉晶界的擴散率和擴散能壘。結果表明，H在(001)表面的擴散比在塊體材料和晶界處快。由DFT和MD方法估計的H擴散能壘高于由實驗得到的沿晶界的能壘。這表明H在Er2O3內的擴散由晶界控制。

2.2 有限元方法及其在腐蝕研究中的應用

有限單元法是隨著電子計算機的發展而迅速發展起來的一種現代計算方法。它是50年代首先在連續體力學領域-飛機結構靜、動態特性分析中應用的一種有效的數值分析方法，隨后很快廣泛的應用于求解熱傳導、電磁場、流體力學等連續性問題。有限元法是一種為求得偏微分方程邊值問題近似解的數值技術。它通過變分方法，使得誤差函數達到最小值并產生穩定解。

ABAQUS是學術研究中使用最為廣泛的有限元分析軟件之一。ABAQUS基于其豐富的單元庫，可用于模擬各種復雜的幾何形狀，并且擁有豐富的材料模型庫，可用于模擬絕大多數常見的工程材料，如金屬、聚合物、復合材料、橡膠、可壓縮的彈性泡沫、鋼筋混凝土以及各種地質材料。隨著實驗室實驗方法得到很大發展的同時，有限元模擬方法也取得了很大的進展，除一般的線性和非線性有限元模擬，還出現一些新的有限元模擬技術。本文主要介紹內聚力單元(Cohesive elements)的有限元模擬方法。

ABAQUS提供一種內聚力單元，用以模擬兩個部分之間的粘性連接，一般來說，它要求粘結材料尺寸和強度都小于粘結部分(比如多層復合材料的膠粘層)，進而可以利用內聚力單元模擬材料的斷裂。

這個單元似乎與C3D8實體單元的結構一致，但最大的區別在于橫縱向尺寸的比例，C3D8單元在模擬這種大橫縱比結構的時候，已經無法給出精確的解答了。這就是內聚力單元存在的意義，此外，從內聚力單元的變形機理來看，中面雖然能夠承受拉伸和剪切應變，但并不能產生任何應力。因此，內聚力單元只能支持垂直于上下表面的牽引-分離破壞準則。當需要模擬諸如扁平結構破壞或者某結構沿一個平面撕裂這種工況的時候，內聚力單元提供了一種無需細化網格的解決方案，可以有效的降低網格數量，提高運算效率。

基于力學的有限元模擬在腐蝕中應用的研究并不是很多，大多數研究集中在材料表面和內部裂紋萌生和擴展問題，尤其是氫敏感體系。Simonovski等[35]提出了利用內聚單元法在晶粒尺度模擬晶間開裂行為的方法。該方法將0厚度的內聚力單元嵌入晶粒之中代表晶界。利用traction-seperation法則進行應力和晶界損傷的計算，結果表明，該方法能夠有效的評價多晶體的晶間損傷和開裂行為。Scheider等[36]用內聚力有限元模型研究了氫擴散對裂紋擴展的影響。Pouillier等[37]利用內聚力模型模擬鋁合金的氫致晶間開裂行為，所得結果與電子背散射技術測得的結果一致。

綜上所述，計算模擬方法在腐蝕領域得到了廣泛的應用，并取得了突破性的進展。計算模擬可為實驗結果提供理論指導，也可為數據分析提供依據。因此，計算模擬方法是研究腐蝕行為與機理的有效手段。

2.3 跨尺度計算與模擬方法研究進展

金屬結構材料在環境中腐蝕的過程從表面對腐蝕性離子的吸附、點蝕與微裂紋的萌生到點蝕坑與裂紋的形成，整個過程從肉眼難以察覺的原子尺度到可觀察的宏觀尺度。通常來講，受到反應時間和尺度的限制，原子尺度的反應很難利用實驗手段進行觀測，這也導致原子尺度的計算與模擬結果很難得到驗證。因此，建立原子尺度到宏觀尺度的跨尺度的計算方法，將原子尺度的輸出結果輸入到宏觀尺度的計算之中，這樣既解決了小尺度實驗難以實現的問題，又可對宏觀尺度的結果進行實驗驗證進而間接驗證原子尺度的計算結果。除此之外，腐蝕環境復雜多變，利用跨尺度的計算與模擬方法能夠對單一因素的影響進行系統研究。

近年來，跨尺度計算已經出現在不同領域的研究之中。多尺度的計算方法一般可分為尺度分離和尺度間耦合兩種，前者著眼于采用不同尺度來分析對象的不同部分，后者注重尋找尺度之間的關聯[38]。通常來講，尺度之間的劃分為宏觀(>10-4m)、介觀(介于10-6m與10-4m之間)、微觀(介于10-7m到10-6m)和納觀(<10-7m)[39]。對于腐蝕計算而言，既需要利用不同尺度的計算方法研究腐蝕發生的不同階段，也需要探究不同尺度之間的耦合關系，通過選擇合適的計算方法對腐蝕的行為與機理進行研究。目前，常見的跨尺度計算模型大多集中在對材料的塑性形變和斷裂過程的研究，如本課題組[40]利用第一性原理與內聚有限元相結合的跨尺度計算方法研究了鋁合金的晶間開裂行為。首先，利用第一性原理研究了H原子在Al內部和晶界處的吸附、溶解和擴散行為，結果表明H原子優先溶解在晶界間隙位且晶界為H原子的快速擴散通道；其次，將不同濃度的H原子引入到晶界間隙位并計算其晶界結合能，結果表明晶界結合能，即晶界強度隨H原子濃度的增加而線性下降；最后，將第一性原理計算獲得的晶界結合能作為斷裂能輸入到內聚有限元計算之中，模擬鋁合金的晶間開裂行為。

3 鋁合金腐蝕與氫致開裂的模擬計算

3.1 鋁鈍化膜點蝕萌生機理

鋁表面覆蓋的氧化物在溶液中或在潮濕空氣中，會吸附環境中的H2O分子形成羥基化的水合表面。當溶液中存在Cl-時，它會與溶液中的OH-在鈍化膜表面發生競爭吸附形成Me-Cl-或是MeO(H)-Cl-的產物，產物的生成導致鈍化膜中的金屬離子與表面的結合力降低從而易于離開表面溶于溶液之中。隨著吸附反應的發生，更多的金屬陽離子遷移到溶液中導致鈍化膜快速減薄，降低其耐蝕性能。關于鈍化膜的局部減薄造成點蝕的機理已被廣泛討論過，但尚未達成一致的說法，尤其是在原子尺度的機理研究還未見報道。本部分利用基于密度泛函的第一性原理計算研究Cl-在金屬鋁鈍化膜表面的吸附和向亞表面滲透的過程。通過對取代和滲透反應的構型、能量和能帶結構的分析來探究Cl-與鈍化膜反應的本質及其導致鈍化膜破裂的機理。

圖1a為Cl-取代羥基化γ-Al2O3(010)表面OH-的構型圖。隨著數值OH-/Cl-的減小，表面發生不同程度的重構。對于OH∶Cl>8∶8，Cl-取代的OH-來自解離的H2O，對氧化膜本身的結構未產生明顯的影響。對于OH∶Cl=2∶14，有14/16的表面OH-吸附位被取代，由構型圖可知表面第一層Al原子均被拉出表面，發生明顯的重構現象。對于OH∶Cl=0∶16來說，表面OH-均被Cl-取代，第一層表面原子與Cl-形成腐蝕產物。圖1b為隨表面Cl-濃度的增加Cl-取代OH-的能量變化，當表面Cl-濃度較低時(小于OH∶Cl=12∶4中所示濃度)，隨著表面Cl-濃度的增大取代能負移，說明Cl-濃度的增加促進其在表面的吸附；當表面Cl-濃度介于OH∶Cl=12∶4和OH∶Cl=8∶8之間時，隨著Cl-濃度的增加取代能正移，說明當表面存在一定吸附Cl-時，后續的Cl-吸附反應難度增加；當表面Cl濃度大于OH∶Cl=8∶8所示，Cl開始取代表面OsH，取代的難度進一步加大。造成能量正移的原因為表面吸附Cl-之間的靜電斥力。當表面全部被Cl-取代時，取代能略有降低，但仍明顯高于低Cl-濃度時。這是因為全部取代使得表面處于對稱結構，降低了靜電斥力間的不均勻，使得吸附結構變得穩定。綜上所述，吸附的Cl-和OH-對后續Cl-的取代反應均有阻礙作用，控制吸附反應發生的因素為表面吸附離子之間的靜電斥力。

Cl-取代表面OH-的反應為放熱反應，隨著取代Cl-濃度的增加，取代能正移，反應驅動力減小。若Cl-取代羥基化γ-Al2O3(010)表面的吸附的來自H2O解離的OH-，則對鈍化膜的結構無明顯影響；若取代鈍化膜本身與H+結合成的OH-時，則發生明顯結構變化，隨著取代Cl-濃度的增大，表面Al與Cl-結合生成無序氯化物層。

圖1 Cl-在水合γ-Al2O3(010)表面發生取代反應：(a)原子構型圖; (b)能量隨表面Cl-濃度的變化Fig.1 Substitution reaction between chloride ions and hydrous γ-Al2O3(010) surface: (a)atomic configuration diagram; (b)variation of energy as chloride ions’ concentration

3.2 氫原子的吸附、溶解和擴散

H原子能夠在鋁表面吸附并進入材料內部，溶解在間隙位。Al (111)表面共有4類對稱吸附位，分別為fcc位、hcp位、top位和bridge位。通過計算表明單一H原子僅能吸附在表面top位和fcc位，當H原子在hcp和bridge位時會自動遷移到最近鄰fcc位。對于H原子在材料內部的溶解，本節共研究了Al(111)層狀表面的5種間隙溶解位，如圖2所示，分別為表面與亞表面之間的兩類四面體間隙位和一類八面體間隙位；塊體材料內部的四面體和八面體間隙位。

圖3a和3b為H原子在晶界∑5(100)/36.87°不同溶解間隙位的遷移路徑和Al原子空位的位置，其中I3GB為靠近晶界的晶內溶解位，I4-1GB和I4-2GB為等效最優溶解位。本節選取了兩條擴散路徑，分別為I1GB→I4-1GB和I5GB→I4-1GB。由圖3c和3d，對于互相垂直的兩條遷移路徑I1GB→I4-1GB和I5GB→I4-1GB，Al原子空位的存在并未降低其擴散能壘，而是大幅度的提高了擴散能壘。當不存在Al原子空位時，I1GB→I4-1GB的遷移能壘為12.96 kJ/mol。Al原子空位的存在使其遷移能壘提高到76.41 kJ/mol，能壘大幅度升高63.45 kJ/mol。對于I5GB→I4-1GB的遷移路徑，無空位存在時的能壘為9.31 kJ/mol，空位的出現使其能壘降低到3.55 kJ/mol，說明Al空位促進H原子沿該路徑的遷移。由以上研究可知，Al原子空位的存在對H原子在晶界的遷移具有很大的影響，Al原子空位能夠促進晶內的H原子遷移進入晶界；Al原子空位對晶界內H原子遷移的作用具有方向性，在某些方向能夠促進H原子的遷移，而在某些方向卻大大抑制H原子的遷移。

圖2 H原子沿Al(111)表面吸附和向內部溶解的間隙位：(a)TIS1-1，表面與亞表面之間的第一類四面體間隙；(b)TIS1-2，表面與亞表面之間的第二類四面體間隙位；(c)OIS1，表面與亞表面之間八面體間隙位；(d)TIS2，金屬Al內部四面體間隙位；(d)OIS2，金屬Al內部八面體間隙位[40]Fig.2 Interstitial sites of H atom adsorbed and dissolved along Al(111) surface: (a) TIS1-1, the first type of tetrahedral interstitial site between surface and subsurface; (b) TIS1-2, the second type of tetrahedral interstitial site between surface and subsurface; (c) OIS1, octahedral interstitial site between surface and subsurface; (d) OIS2, octahedral interstitial site in Al bulk[40]

3.3 氫致鋁晶間開裂行為的跨尺度計算

利用Neper軟件產生多晶體結構模型，如圖4所示。晶粒之間的界面利用FORTRAN語言編寫的程序產生。多晶體模型包含300個晶粒，1387個晶界面。線性四面體單元(ABAQUS 類型 C3D4) 被用于劃分晶內部分網格, 線性三角單元(ABAQUS類型 COH3D6) 被用于劃分晶界部分網格。 多晶體晶粒的取向為正交各向異性。沿x方向的10-3mm的位移邊界條件被施加在右側表面，左側表面的節點沿x、y、z三個方向固定。

圖3 H原子在晶界(100)/36.87°，∑5的溶解位示意圖：(a)側視圖；(b)俯視圖(其中箭頭表示為不同的遷移路徑，黃色四邊形標記處為Al原子空位的位置)。H原子在鋁晶界間隙位擴散的路徑：(c)原子在有無空位情況下沿I1GB→I4-1GB的遷移路徑；(d)H原子在有無空位情況下沿I5GB→I4-1GB的遷移路徑(綠色箭頭表示的數值為有無空位情況下遷移能壘的差值)Fig.3 Dissolved sites of H atom in (100)/36.87°，∑5 interface: (a) side view; (b) top view; (c) migration path of H atom along I1GB→I4-1GB with and without vacancies; (d) migration path of H atom along I5GB→I4-1GB with and without vacancies. Green arrows represent the difference values between migration energy barrier with vacancies and without vacancies

圖4 多晶體模型圖Fig.4 Polycrystal model diagram

利用第一性計算原理方法研究晶界偏聚對晶界強度的影響，將不同濃度H原子置于Al∑5 [100]/36.87°晶界超胞穩定吸附位，計算其晶界結合能，表明晶界結合能隨H濃度升高而線性降低，如圖5所示，即晶界強度線性降低。H原子對晶界的弱化效應歸因于原子H的存在弱化了晶界處Al原子之間的金屬鍵。

圖5 Al∑5 [100]/36.87°晶界結合能隨H濃度增加的變化Fig.5 The relationship between H concentration and cohesive energy of Al GB ∑5 [100]/36.87°

當一定量的H原子從環境中進入材料內部，H原子能夠偏聚在晶界，導致原子之間結合力降低，這是鋁合金失效的主要原因之一；雖然局部晶界應變已能夠通過實驗測試獲得，但應力的分布仍然難以測得,內聚力有限元方法能夠對晶界的應力應變分布和損傷行為進行細致的計算與模擬，其中最重要的輸入參數—晶界結合能，能夠通過第一性原理計算獲得。通過內聚有限單元法對鋁多晶體模型的晶間開裂行為進行模擬，其中利用第一性原理計算得到的晶界結合能將作為斷裂能這一參數輸入到內聚有限元計算中，實現原子尺度與宏觀尺度之間的銜接，從而實現跨尺度計算。圖6為晶界單元的損傷云圖，可以看出，隨著晶界結合能的降低，失效單元增多，隨著失效單元的累積使得微觀晶間出現微觀裂紋。

圖6 不同晶界結合能下內聚有限單元法計算得到的SDEG損傷分布云圖[40]Fig.6 SDEG damage distribution diagrams of polycrystal model with different binding energy of grain boundary via cohesion finite element method[40]

4 結 語

本文對計算模擬方法在金屬腐蝕研究中的應用與進展進行了綜述，將第一性原理與有限元計算相結合，建立了從原子尺度到宏觀尺度的跨尺度計算方法。本文以鋁及鋁合金為例，利用第一性原理計算了Cl與表面鈍化膜的相互作用來研究鈍化膜的點蝕機理。本文還建立了Al-H體系，研究了H原子在Al表面和內部的吸附、溶解，擴散與偏聚，并計算H原子偏聚對晶界強度的影響。最后將原子尺度第一性原理計算獲得的能量輸入到宏觀尺度的有限元計算之中模擬開裂行為。本研究基本完成鋁大氣腐蝕行為與機理的跨尺度計算與模擬框架的搭建，此方法適用于服役于大氣環境中具有Cl-點蝕敏感性和氫脆敏感性的結構材料，對工程應用具有指導作用。

本課題組在后續的工作中還針對金屬鈍化膜的性能開展了大量的研究，利用第一性原理與電化學測試相結合的方法研究鈍化膜的穩定性[41,42]及其半導體性質對其耐蝕性能的影響。

References

[1] Xiao Jimei(肖紀美), Cao Chunan(曹楚南).PrinciplesofMaterialCorrosion(材料腐蝕學原理)[M]. Beijing: Chemical Industry Press, 2002.

[2] Zhu X, Zi G.Construction&BuildingMaterials[J], 2017, 137: 330-344.

[3] Maulana A, Dipoyono H K, Khairurrijal K.IndonesianJournalofPhysics[J], 2015, 18(1): 5-9.

[4] Wang H, Han E H.ElectrochimicaActa[J], 2013, 90(5): 128-134.

[5] Laycock N J, Krouse D P, Hendy S C,etal.ElectrochemicalSocietyInterface[J], 2014, 23(4): 65-71.

[6] Wang H T, Han E H.Materials&Corrosion[J], 2015, 66(9): 925-930.

[7] Harrison T J, Crawford B R, Brandt M,etal.ComputationalMaterialsScience[J], 2014, 84(1): 74-82.

[8] Li Q K, Zhang Y, Shi S Q,etal.MaterialsLetters[J], 2012, 56(6): 927-932.

[10] Zhang G A, Zeng L, Huang H L,etal.CorrosionScience[J], 2013, 77: 334-341.

[11] Liu Zhixiao(劉智驍).TransactionsofNonferrousMetalsSocietyofChina(中國有色金屬學報)[J], 2013, 23(4): 1160-1167.

[12] Pan Changchang(潘昌昌). Dissertation for Master[D]. Lanzhou：University of Lanzhou Technology, 2016.

[13] Venkatraman M S, Cole I S, Emmanuel B.ElectrochimicaActa[J], 2011, 56(20): 7171-7179.

[14] Indira K, Nishimura T.JournalofMaterialsEngineering&Performance[J], 2016, 25(10): 1-14.

[15] Giner I, Keller A.CorrosionScience[J], 2015, 100: 496-503.

[16] Ceylan S, Sergiy B, Martin S,etal.CorrosionScience[J], 2012, 58(5): 307-314.

[17] Di Stefano D, Mrovec M, Els?sser C.ActaMaterialia[J], 2015, 98: 306-312.

[18] Kulkova S E, Kulkov S S, Bakulin A V,etal.InternationalJournalofHydrogenEnergy[J], 2012, 37(8): 6666-6673.

[19] Hohenberg P, Kohn W.PhysicalReview[J], 1964, 136：B864-B871.

[20] Kohn W, Sham L J.PhysicalReview[J], 1965, 140(4A): A1133-A1138.

[21] Perdew J P, Burke K, Ernzerhof M.PhysicalReviewLetters[J], 1996, 77(18): 3865-3868.

[22] Wei Xin, Dong Chaofang, Chen Zhanghuaetal.TransactionsofNonferrousMetalsSocietyofChina[J], 2015, 25: 4102-4109.

[23] Guo F Y, Long C G, Zhang J,etal.AppliedSurfaceScience[J], 2015, 324: 584-589.

[24] Carrasco J, Hodgson A, Michaelides A.NatureMaterials[J], 2012, 11(8): 667-674.

[26] Tonigold K, Gross A.JournalofComputationalChemistry[J], 2012, 33(6): 695-701.

[27] Wei X, Dong C F, Chen Z H,etal.RSCAdvances[J], 2016, 6:79836-79843.

[28] Wei X, Dong C F, Chen Z H,etal.RSCAdvances[J], 2016, 6: 56303-56312.

[29] Digne M, Sautet P, Raybaud P,etal.JournalofCatalysis[J], 2012, 226(1): 54-68.

[30] Wei X, Dong C, Chen Z,etal.AppliedSurfaceScience[J], 2015, 355(15): 429-435.

[31] Fernandez N, Ferro Y, Kato D.ActaMaterialia[J], 2015, 94(1): 307-318.

[32] Kristinsdóttir L, Skúlason E.SurfaceScience[J], 2012, 606:1400-1404.

[33] Ferrin P, Kandoi S, Nilekar A U,etal.SurfaceScience, 2012, 606(7-8): 679-689.

[34] Mao W, Chikada T, Suzuki A,etal.JournalofNuclearMaterials[J], 2014, 455(1): 360-365.

[35] Simonovski I, Cizelj L.EngineeringFractureMechanics[J], 2013, 110: 364-377.

[36] Zhu L K, Yan Y, Qiao L J,etal.CorrosionScience[J], 2013, 77(77): 360-368.

[37] Pouillier E, Gourgues A F, Tanguy D,etal.InternationalJournalofPlasticity[J], 2012, 34(34): 139-153.

[38] Lu Xinzheng(陸新征), Lin Xuchuan(林旭川), Ye Lieping(葉列平).JournalofHuazhongUniversityofScienceandTechnology(華中科技大學學報)[J], 2008, 25(4): 76-79.

[39] Guo Yafang(郭雅芳), Wang Chongyu(王崇愚).MaterialsReview(材料導報)[J]，2001，15(7)：9-11.

[40] Wei X, Dong C F, Chen Z H,etal.RSCAdvances[J], 2016, 6:27282-27292.

[41] Xu A N, Dong C F, Wei X,etal.ElectrochemistryCommunications[J], 2016, 68: 62-66.

[42] Kong D C, Xu A N, Dong C F,etal.CorrosionScience[J], 2017, 116: 34-43.