許競翔，馮建方，孫 強，褚振華

(上海海洋大學工程學院，上海 201306)

0 前 言

海洋環(huán)境十分惡劣，金屬表面水含量、溫度、輻照度、氯離子濃度、SO2等污染物等環(huán)境因素對金屬材料的腐蝕有很強的影響，從而對其力學性能和幾何形狀等造成不同程度的影響[1，2]，造成設(shè)施結(jié)構(gòu)損傷，使役壽命縮短，嚴重時危害生產(chǎn)安全，并造成相當大的財政負擔[3]，其中船舶和海洋工程裝備中金屬材料的腐蝕防護成本占比很高，對其腐蝕機理的研究尤顯重要。

長期以來，有關(guān)金屬材料腐蝕的研究大多數(shù)是通過試件實驗進行的，主要是通過模擬實際環(huán)境條件，觀察金屬材料的特征和變化，從腐蝕產(chǎn)物、腐蝕表面狀態(tài)、腐蝕產(chǎn)物的演變等方面來研究各種金屬材料在不同環(huán)境中的腐蝕機制和原理[4-8]。但是金屬發(fā)生腐蝕的環(huán)境十分復雜，受到眾多因素以及因素之間耦合作用的影響，任何單一因素對腐蝕行為的影響都可能受其他因素的干涉，腐蝕的機理不太清晰，且腐蝕實驗所需成本往往較高。數(shù)值計算模擬相對成本較低，該方法從原子/分子尺度研究腐蝕在微觀尺度上的動態(tài)變化過程，如金屬與環(huán)境的化學反應、化學鍵的斷裂和形成、表面的缺陷和擴散、應力應變的變化等等。

隨著計算機技術(shù)的快速發(fā)展，數(shù)值模擬已經(jīng)成為研究腐蝕等領(lǐng)域的方法之一，其中分子動力學和密度泛函(DFT)理論等是主要的數(shù)值模擬方法[9-13]。金屬的腐蝕是指金屬和環(huán)境之間發(fā)生化學或電化學反應，金屬成為離子態(tài)或與其他粒子結(jié)合形成化合物，其實質(zhì)是原子間化學鍵的斷裂和形成，現(xiàn)有的研究[9-13]充分證明分子動力學有能力描述金屬的腐蝕過程。分子動力學模擬是通過構(gòu)造出比較理想的模型，改變單一變量(如溫度、表面取向、氯離子濃度等)模擬在真實環(huán)境中金屬材料腐蝕所發(fā)生的動態(tài)過程，然后與實驗研究對比，明確腐蝕機理，從而找到更加高效的抑制或減緩腐蝕發(fā)生的方法[13]。分子動力學模擬已成為研究金屬材料腐蝕的重要方法，目前已成功用于摩擦磨損[14]、超精密加工[15]、微納米元器件[16]等研究領(lǐng)域。

1 原子間的勢模型

分子動力學算法的基礎(chǔ)，在于從原子間相互作用勢函數(shù)出發(fā)，基于牛頓第二定律來求解體系中各原子的運動軌跡，進而獲取體系隨時間的演化過程。原子間相互作用勢函數(shù)的可靠性直接決定所得原子運動軌跡的準確性，從而在分子動力學模擬中具有非常核心的作用[17-20]。

在多粒子系統(tǒng)中，牛頓第二運動定律的二階微分方程定義為[21，22]：

Fi=miai

(1)

式中：Fi是原子i的力,mi和ai分別是原子i的質(zhì)量和加速度。

用有限差分法對式(1)中的各自由度進行時間積分：

Fi=-?iV

(2)

式中：V是原子間勢能函數(shù)，F(xiàn)i是其梯度。

原子間勢能函數(shù)的計算是分子動力學模擬中計算量最大的一步，通常占整個模擬時間的95%[21]。描述原子間相互作用的勢函數(shù)很多，早期主要是通過兩體勢函數(shù)[23]描述，如Lennard - Tones(LJ)勢和Morse勢等。但是兩體勢函數(shù)在描述原子間化學鍵的斷裂和形成方面還有些不足。隨著研究和發(fā)展，出現(xiàn)了考慮角度效應和靜電勢等的多體勢函數(shù)，而描述化學反應較好的勢函數(shù)主要有修正的嵌入原子法[23，24](Modified Embedded Atom Method，MEAM)和反應力場(Reaction Force Field，ReaxFF)勢函數(shù)等等[25]。

1.1 MEAM勢函數(shù)

嵌入原子法(Embedded Atom Method，EAM)勢函數(shù)主要描述的是所有原子核之間的排斥作用和原子的原子核與核外電子及原子周邊其他電子的靜電作用，具體形式通過擬合完整純金屬體系的點陣常數(shù)、升華能、彈性常數(shù)、空位形成能、合金潛熱等基本物理量來定義[23]。EAM勢函數(shù)定義為[24]：

(3)

式中：EEAM表示系統(tǒng)的總能量，i和j分別表示兩原子，E表示將原子i嵌入到電子密度為ρi的位置處i所需要的嵌入能，φ表示相距r的i和j原子之間的兩體勢。

而在修正的嵌入原子法(MEAM)中，系統(tǒng)的總能量為嵌入原子勢和靜電勢(Electrostatic Potential，ES)之和。MEAM勢函數(shù)定義為[24]：

Esystem=EEAM+EES

(4)

式中：Esystem表示修正后的系統(tǒng)的總能量，EEAM和EES分別表示嵌入原子勢能和靜電勢能。EES定義為[24]：

(5)

式中：Zj表示有效的原子核電荷，一階導數(shù)χ是電負性，二階導數(shù)J與自庫侖斥力相關(guān)，ωi(rij)表示短距離核吸引力積分，μij(rij)表示短距離庫侖相互作用積分，qi是離子i的電荷。

Lee等[26]通過MEAM開發(fā)了Fe - Al二元體系的原子間勢，該電勢可以描述Fe - Al二元合金的各種基本物理特性 - 結(jié)構(gòu)、彈性和熱力學特性、缺陷形成行為以及缺陷之間的相互作用，且這些性能與實驗數(shù)據(jù)和第一性原理計算的數(shù)據(jù)基本吻合，證明了該方法從原子角度研究各種缺陷形成行為及其對高鋁鋼以及Fe - Al二元合金的力學性能的影響方面的適用性。

Jelinek等[27]采用MEAM開發(fā)了鋁、硅、鎂、銅和鐵的相互作用勢，根據(jù)從頭算(ab initio method，基于量子力學基本原理直接求解薛定諤方程的量子化學計算方法)得到的NaCl參考結(jié)構(gòu)中元素對的結(jié)構(gòu)和彈性特性，建立了元素對的MEAM參數(shù)，并進行了調(diào)整，重現(xiàn)了最穩(wěn)定二元化合物的從頭計算生成熱。然后對幾種二元化合物的缺陷形成能、平衡體積、彈性模量和生成熱與從頭算模擬和實驗結(jié)果比較，驗證了MEAM勢的可行性。

Campbell等[24]通過MEAM模擬半徑為10 nm的鋁納米團簇的氧化過程，基于電負性均衡原理考慮了鋁和氧之間電荷轉(zhuǎn)移的影響。鋁向表面擴散，而氧向團簇內(nèi)部擴散，鋁的擴散率比氧在氧化物中的擴散率高30%～60%。氧化滲濾的過程分為3步，分別為從氧氣分子吸附到鋁表面、氧原子通過fcc鋁的四面體和八面體的位置向基體擴散、最后形成[OAl4]四面體簇，氧化膜厚度隨時間的變化如圖1所示，該氧化過程解釋了在典型模擬的早期階段氧化物生長的減速，同時也解釋了為何氧化層的厚度和平均溫度隨時間的延長而不斷增加。

1.2 ReaxFF勢函數(shù)

反應力場(ReaxFF)是專門設(shè)計出來描述化學反應性、化學鍵的斷裂和形成、表面、缺陷、擴散等的一種特殊類型的力場[25]，它同時考慮了非鍵合和鍵合的能量函數(shù)。它通過能量描述化學鍵，有著較高的計算效率，彌補了量子化學方法和傳統(tǒng)分子動力學之間的差距[25]。

ReaxFF原子間勢能函數(shù)定義為[25，28，29]：

Esystem=Ebond+Eover+Eunder+Eval+Etor+Elp+
EH-bond+Evdwaals+Ecoulomb

(6)

式中：Esystem表示系統(tǒng)的總能量，Ebond、Eover、Eunder、Eval、Etor、Elp、EH-bond、Evdwaals和Ecoulomb等能量項分別對應鍵能、欠配位的能量矯正項、過配位的能量矯正項、價角能、扭轉(zhuǎn)角能、孤對能、氫鍵能、范德華能和庫侖能。最后兩個能量項范德華能Evdwaals和庫侖能Ecoulomb計算所有原子對，而其他所有能項均由計算粒子i和j之間的鍵級獲得。

原子電荷是使用與幾何有關(guān)的電荷分布來計算的，該分布通過電荷平衡法確定[30-32]。單個原子電荷和靜電能E(q)隨時間動態(tài)地變化。相應的靜電能E(q)定義為：

(7)

式中：qi是離子i的電荷；N是離子的總數(shù)；q是電荷長度N的矢量；wi和ui分別是離子i的電負性和硬度。Tap(rij)是一個七階錐度函數(shù)，當帶電粒在非鍵合截止半徑內(nèi)或外移動時，避免能量不連續(xù)；Kc是介電常數(shù)；rij是原子i和j之間的距離。原子電荷與系統(tǒng)的幾何位置有關(guān)，在每一時間步上都進行計算。

Nabankur等[33]利用ReaxFF分子動力學模擬研究了陽離子Li+、N a+、K+和陰離子Cl-對水的結(jié)構(gòu)和動力學特性的影響，且已經(jīng)成功地描述了這些系統(tǒng)固有的化學反應和擴散行為，能夠?qū)滏I動力學的定性趨勢與不同體系的物種分布聯(lián)系起來，使得其他電解質(zhì) - 水系統(tǒng)和熱力學狀態(tài)的模擬成為可能。

Gabriele等[34]通過ReaxFF描述了銀 - 氧和銀 - 硫反應機理，Ag - O和Ag - S的反應快照如圖2所示。

圖2中下方基底處較大的原子代表銀原子，圖2a中中上方小原子代表硫原子，圖2b中中上方小原子代表氧原子。由圖2可知Ag - O和Ag - S具有不同的反應機理和反應速率。Ag2O的生長必須通過氧化物表面的O2解離和O通過氧化物層向氧化物/Ag間相擴散進行，O2的解離黏著系數(shù)隨著氧覆蓋率的增加而急劇下降。一旦S8分子接近銀，二者就不可避免地結(jié)合并發(fā)生反應，而且在第一層硫化銀形成之后，硫化銀的形成也相對較快，這是因為硫化物在銀表面以外的生長不需要缺陷的存在。這個模擬證明了銀極易受到硫腐蝕(失去光澤)，而很難與O2發(fā)生反應。

2 金屬腐蝕分子動力學模擬模型

建立合適的金屬腐蝕分子動力學反應模型是得到正確和可靠的模擬結(jié)構(gòu)的首要條件。目前用于金屬腐蝕的分子動力學模擬的反應模型主要有3種，分別為金屬團簇 - 溶液/氣體反應模型、溶液/氣體 - 金屬表面 - 溶液/氣體反應模型和金屬表面 - 溶液/氣體反應模型。

2.1 金屬團簇 - 溶液/氣體反應模型

圖3為金屬團簇 - 溶液/氣體反應模型的截面示意圖。在該模型中，溶液/氣體中的原子與金屬團簇表面直接接觸反應。Li等[35]使用金屬團簇 - 溶液/氣體反應模型研究鋁團簇與水分子的反應，采用分子動力學的方法模擬了不同數(shù)量的水分子與鋁團簇的反應行為。

2.2 溶液/氣體 - 金屬表面 - 溶液/氣體反應模型

圖4為溶液/氣體 - 金屬表面 - 溶液/氣體反應模型的示意圖。在該模型中，溶液/氣體中的原子與金屬不同表面取向的表面(雙面)直接接觸反應。Liu等[36]使用溶液/氣體 - 金屬表面 - 溶液/氣體模型研究鋁(100)與O2的反應，采用分子動力學的方法模擬了鋁的氧化過程。

2.3 金屬表面 - 溶液/氣體反應模型

圖5為金屬表面 - 溶液/氣體反應模型的示意圖。在該模型中，溶液/氣體中的原子與金屬不同表面取向的表面(單面)直接接觸反應，如圖5所示，在圖5的左邊固定一部分金屬基體，在圖5的右邊一定距離處設(shè)置一個反射墻，在兩端分別預留一定距離的真空。這種設(shè)置的目的是讓Z方向展現(xiàn)為非周期性。DorMohammadi等[37]使用金屬表面 - 溶液/氣體反應模型研究了鐵基體不同的表面取向的表面與純水的反應，采用分子動力學的方法模擬了鐵基體在純水中腐蝕的初始階段。

從以上分析可以看出，對于不同的金屬材料體系，需要按要求采取不同的反應模型。一般來說，在研究金屬材料腐蝕的過程中，大多數(shù)采用金屬表面 - 溶液/氣體反應模型[25-27，32，34]，因為在實際中，金屬表面往往只有一層與環(huán)境(溶液)直接接觸，所以金屬表面 - 溶液/氣體反應模型被認為最能真實反映宏觀尺度的金屬材料腐蝕。但是目前如何橋接原子尺度與真實宏觀尺度成為一個重要的問題[35]，這也是未來金屬腐蝕分子動力學研究的重點方向之一。

3 金屬腐蝕分子動力學模擬研究進展

海洋工程裝備常用金屬材料有鋼、鐵、鈦合金、銅合金、鋁合金、鎂合金等[38]。海洋中的海水是一種具有一定流速的腐蝕性非常強的天然電解質(zhì)溶液[38]。國內(nèi)外研究表明[23，38-40]，金屬表面水含量、溫度、氯離子濃度、表面取向、外電場等因素對海洋金屬材料腐蝕都有著顯著的影響，且這些因素之間有耦合作用，使得腐蝕變得復雜化。但是分子動力學模擬可以確定單一因素對金屬材料腐蝕的影響。

3.1 金屬表面水含量

一般將海洋腐蝕區(qū)域分為海洋大氣區(qū)、浪花飛濺區(qū)、潮差區(qū)、海水全浸區(qū)和海底泥沙區(qū)[38]，在這些區(qū)域中金屬材料表面有一層相對較厚的水膜。當金屬材料表面的相對濕度達到臨界相對濕度時，腐蝕速度將大幅上升[40]。目前采用分子動力學模擬的方法對金屬表面水含量對金屬材料腐蝕的影響的報道已經(jīng)有很多，但是尚沒有形成統(tǒng)一的結(jié)論。

一些分子動力學模擬研究表明，部分金屬材料如鐵、鋁等在純水中會被腐蝕。DorMohammadi等[37]分別建立了鐵(100)、(110)、(111)與純水的金屬溶液模型，圖6是鐵(110)基體與純水的初始構(gòu)型，圖6中左邊整齊排列的是鐵原子，右邊是水分子團簇。

模擬結(jié)果表明鐵基體在初始階段會被純水腐蝕，且鐵(100)、(110)、(111)基體表面隨著反應的進行將生成鐵的不同價態(tài)的氧化物，并有鐵原子溶解到純水中。Li等[35]建立了鋁納米團簇與純水的金屬團簇 - 溶液/氣體反應模型，模擬結(jié)果表明鋁團簇與純水發(fā)生反應，鋁的表面會生成氧化鋁并伴隨著氫氣生成，在一段時間后，鋁表面氧化膜覆蓋團簇表面，鋁團簇不再與水反應。

另一些分子動力學模擬研究表明，鎳等金屬材料不會被水腐蝕。Assowe等[25，29]建立了鎳(111)與純水的金屬溶液模型，模擬結(jié)果表明，在沒有外電場的情況下，鎳(111)基體表面不和水反應，也沒有鎳原子溶解到水中。

從以上分析可以看出，金屬表面水含量會影響部分金屬的腐蝕。一些本身不與水反應的金屬，在加入了外電場等其他的影響因素后被腐蝕。

3.2 溫 度

溫度升高，原子吸收能量，相對而言更活躍，原子之間化學反應速率變快，因此溫度對金屬材料的腐蝕也會產(chǎn)生顯著的作用。地球表面的太陽輻射能大部分被海洋所吸收，晝夜交替、洋流、大氣環(huán)流等使得海水的溫差很大。在海洋工程裝備工作的時候，局部溫度可能很高。而受限于實驗設(shè)備，海洋工程材料在極端溫度下的腐蝕缺乏系統(tǒng)研究。分子動力學模擬可以很穩(wěn)定地控制系統(tǒng)中的溫度變化，所以分子動力學方法能夠很好地研究溫度對金屬材料腐蝕的影響。

一些分子動力學模擬研究表明，溫度越高，海洋金屬材料腐蝕的速率就越快。DorMohammadi等[37]建立了Fe(110)與純水的金屬溶液模型(如圖6所示)。研究溫度變化范圍為300～350 K，在25 MeV/cm外電場下鐵和水側(cè)氧原子數(shù)的比較，水側(cè)的鐵原子數(shù)隨時間的變化如圖7所示。

從圖7中可以看出，隨著溫度的升高，鐵側(cè)的氧原子數(shù)略有增加，但水側(cè)的鐵原子數(shù)幾乎相同。圖8為不同溫度下150 ps 時25 MeV/cm外電場界面的鐵側(cè)氧原子密度和水側(cè)鐵原子密度與界面Z的關(guān)系圖，結(jié)合圖7a可知，隨著溫度的升高，鐵的氧化略有增加。綜上可得，溫度的升高加快了鐵的腐蝕速率，但是并沒有增加鐵的腐蝕程度。

3.3 氯離子濃度

海洋中含有大量的鹽類，通常把海水近似地看作3.0%(質(zhì)量分數(shù)，下同)或者3.5%的NaCl溶液[39]。但是在實際海洋環(huán)境中，氯離子的濃度是有變化的，在海洋金屬材料表面有一層水膜，水膜的蒸發(fā)和凝聚會導致材料表面鹽分濃縮，改變氯離子的濃度[38]。

一些分子動力學研究表明，金屬的腐蝕和氯離子的濃度有關(guān)。Jeon等[40]建立了 Cu (111)基體在氯離子溶液中的腐蝕模型，圖9為氯離子濃度為10 mol/L的水溶液條件下銅腐蝕模擬的初始模型，研究了水相條件下銅基體與氯的相互作用。300 K下腐蝕模擬弛豫500 ps時，不同氯離子濃度的水性介質(zhì)和銅基底之間的界面的側(cè)視圖如圖10所示，圖10中下方排列規(guī)則的原子代表銅原子，上方緊密排列的有3個原子和2個原子的分子分別代表水分子和氯氣分子。由圖10可知，腐蝕模擬弛豫500 ps時，在0 mol/L和1 mol/L氯離子的模型中，銅離子沒有擴散到溶液中，在10 mol/L氯離子的模型中，銅發(fā)生腐蝕且有銅離子擴散到溶液中。

DorMohammadi等[41]建立了Fe(110)基體在pH=13.5的堿性溶液中的模型，結(jié)果表明在Fe(110)表面形成了一層由富含F(xiàn)e2+的內(nèi)氧化物層和富含F(xiàn)e3+的外氧化物層組成的氧化膜；在堿性溶液中加入氯離子以后，F(xiàn)e(110)表面的氧化膜結(jié)構(gòu)被破壞，圖11為氯化物和氫氧化物與氧化鐵表面相互作用的不同階段。圖11中被方框、菱形、五邊形圈住的小球分別表示鐵原子、氯原子和OH-。如圖11所示，鐵原子以Fe(OH)Cl2和FeCl3的形式與鈍化膜分離，溶解到溶液中；氯化物在去鈍化過程中有2個主要功能：首先，引起薄膜/電解質(zhì)界面附近電解質(zhì)的局部酸化。其次，引起暴露鐵原子的溶解。2種作用都需要電解質(zhì)中存在足夠數(shù)量的氯離子；在電解質(zhì)中使用1 mol/L和2 mol/L NaCl進行的模擬顯示沒有局部酸化和鐵的溶解，表明太低的氯化物濃度不能引發(fā)鐵表面氧化膜的腐蝕。

3.4 表面取向

腐蝕現(xiàn)象的本質(zhì)是金屬材料和環(huán)境中其他物質(zhì)之間的化學反應，即原子間的鍵的斷裂和形成，因此基體的晶體結(jié)構(gòu)、表面原子結(jié)構(gòu)的取向，以及溶液與晶體之間的取向均會對腐蝕產(chǎn)生影響。分子動力學模擬方法在研究表面取向?qū)饘俑g的影響方面有著巨大的優(yōu)勢，國內(nèi)外許多學者對此開展了研究。

一些分子動力學模擬研究表明，在相同條件下，金屬晶體不同的表面取向?qū)饘倩w的特性有不同的影響。Jeon等[42]分別建立了鐵基體不同表面(100)、(110)、(111)的氧化模型，圖12為鐵(100)、(110)、(111)在氧化1 ns時的側(cè)視圖，圖12中下方整齊排列的為鐵原子，雜亂排列的為氧原子。模擬結(jié)果表明，在300 K和900 K下，F(xiàn)e(110)顯示出最多量的氧化物生長，F(xiàn)e(111)次之，F(xiàn)e(100)最少量(如圖12所示)。DorMohammadi等[37]分別建立了Fe(100)、(110)、(111)與純水的金屬溶液模型(如圖6所示)，模擬結(jié)果表明，不同的鐵表面取向?qū)τ阼F的氧化過程和電荷分布沒有明顯的影響。綜合分析得到，鐵基體的不同表面取向不會對鐵的氧化產(chǎn)生顯著的影響。

3.5 外電場

電場是電荷及變化磁場周圍空間里存在的一種特殊物質(zhì)，電場由電荷和隨時間變化的磁場產(chǎn)生，其會對處在場內(nèi)的其他電荷產(chǎn)生力的作用。在微觀尺度下，電場本質(zhì)上是原子內(nèi)原子核和電子之間的吸引力以及原子間的電子相互作用而形成的作用力[43]。反應力場分子動力學模擬在金屬溶液模型之間引入一個外電場，在電場作用下，金屬氧化膜的結(jié)構(gòu)和形態(tài)會發(fā)生變化。

一些分子動力學模擬研究表明，外電場會使本身不與水發(fā)生化學反應的金屬發(fā)生腐蝕。Assowe等[25]建立了Ni(111)與純水的金屬溶液模型，鎳的氧化膜厚度隨電場強度的變化如圖13所示。模擬結(jié)果表明隨著電場強度的增加，氧化膜厚度呈線性增加趨勢，鎳的腐蝕速率加快。

4 結(jié) 語

(1)國內(nèi)外學者對金屬材料腐蝕的反應力場分子動力學模擬進行了很多研究，從這些模擬結(jié)果來看，金屬材料的腐蝕受到多種因素的影響，在分子/原子尺度上的腐蝕和宏觀尺度上的腐蝕存在著明顯的區(qū)別，使得到目前為止仍很難從本質(zhì)上給出腐蝕的起源并歸納普適的腐蝕規(guī)律。

(2)在未來腐蝕分子動力學模擬領(lǐng)域研究中，需要在以下幾個方向開展重點研究：第一，加強對一些不完整金屬氧化膜表面腐蝕機理的研究，例如金屬的點缺陷、位錯、晶界、雜質(zhì)等缺陷模型[44]的研究，從原子角度闡述其腐蝕的機理。第二，建立能夠反應真實宏觀尺寸表面的金屬腐蝕分子動力學模擬的模型，確定在什么參數(shù)下分子/原子尺度可以收斂到真實宏觀尺度。第三，加強金屬材料腐蝕分子動力學模擬結(jié)果與試驗結(jié)果之間的相互關(guān)聯(lián)和橋接，完善腐蝕理論[25]。