金屬表面點(diǎn)蝕的穩(wěn)態(tài)過程數(shù)值研究

崔靜，楊帆，楊霆浩，楊廣峰

中國民航大學(xué)機(jī)場學(xué)院，天津300300

0 引 言

現(xiàn)代艦船裝備結(jié)構(gòu)的系統(tǒng)化和復(fù)雜化程度越來越高，在艦船上使用的新材料也越來越廣泛，因此艦船的腐蝕防護(hù)工作面臨著巨大挑戰(zhàn)。金屬腐蝕會(huì)導(dǎo)致艦船裝備的損壞，造成經(jīng)濟(jì)損失以及腐蝕污染。預(yù)測腐蝕的損傷過程，對于保證艦船裝備的安全和正常運(yùn)行起著重要的作用。艦船所處的海洋環(huán)境復(fù)雜，運(yùn)用實(shí)際測量和經(jīng)驗(yàn)估算等方法對艦船材料表面腐蝕進(jìn)行檢測時(shí)，評估標(biāo)準(zhǔn)不統(tǒng)一，不能實(shí)時(shí)觀測腐蝕發(fā)展過程，而采用計(jì)算機(jī)技術(shù)進(jìn)行數(shù)值模擬，在一定程度上可彌補(bǔ)傳統(tǒng)方法的不足。

腐蝕過程是一個(gè)包含了多組分之間的質(zhì)量傳輸、化學(xué)反應(yīng)或電化學(xué)反應(yīng)、界面動(dòng)力學(xué)等多重因素影響的復(fù)雜過程，而格子Boltzmann方法（Lattice Boltzmann Method，LBM）作為一種介觀方法（Mesoscopic technique），近年來發(fā)展迅速，成功地應(yīng)用于多種微觀相互作用的流體系統(tǒng)問題中。LBM既能兼顧實(shí)際的計(jì)算成本，又能保證模型或方法的合理性。國內(nèi)外運(yùn)用LBM進(jìn)行腐蝕模擬的研究較少，更多學(xué)者采用的是元胞自動(dòng)機(jī)模型（Cellular Automata Model，CAM）進(jìn)行腐蝕模擬，采用CAM進(jìn)行腐蝕模擬對于多組分之間的質(zhì)量傳輸和腐蝕化學(xué)反應(yīng)只進(jìn)行了簡單的元胞轉(zhuǎn)化，并沒有理論支撐。He等［1］首次將化學(xué)反應(yīng)與LBM結(jié)合，解決了流體—壁面表面反應(yīng)與擴(kuò)散問題的耦合。Kang等［2-3］對溶解沉淀問題中單相流體和多相流體的傳輸及流動(dòng)進(jìn)行了研究，將LBM與結(jié)點(diǎn)體積法（Volumn of Pixel，VOP）結(jié)合來模擬物質(zhì)的傳輸過程。Chen等［4］為了解決多相流動(dòng)中的相變和溶解沉淀問題，將LBM和VOP方法在孔隙尺度下相結(jié)合，解決了多相流中的相分離、化學(xué)反應(yīng)、質(zhì)量傳輸和溶解沉淀等問題。張婷等［5］使用MRT-LBM來模擬腐蝕溶液的流動(dòng)，并分析了不同參數(shù)對腐蝕的影響。Liu等［6］使用D3Q19模型模擬了巖石的腐蝕情況。Zhang等［7］在腐蝕模型中加入了電場并與流場成功耦合，更好地模擬了電化學(xué)腐蝕。但是，現(xiàn)有的格子Boltzmann（LB）腐蝕模型中，并未對具體的金屬材料進(jìn)行模擬，更多是使用化學(xué)腐蝕的方法進(jìn)行研究，僅考慮了陽極固體腐蝕，很少涉及陰極反應(yīng)的發(fā)生。

本文將針對現(xiàn)有LB腐蝕模型的不足，采用陰陽極共同作用的方式研究電化學(xué)腐蝕，將使用改進(jìn)后的模型對金屬表面點(diǎn)蝕現(xiàn)象的穩(wěn)態(tài)過程進(jìn)行模擬，得到點(diǎn)蝕坑的形貌變化特征以及各組分濃度的變化情況。通過數(shù)值模擬，分析腐蝕化學(xué)反應(yīng)速率、腐蝕溶液擴(kuò)散系數(shù)、腐蝕產(chǎn)物擴(kuò)散系數(shù)對腐蝕程度的影響。

1 格子Boltzmann腐蝕模型

腐蝕過程包含多相多組分的流動(dòng)與擴(kuò)散、電化學(xué)反應(yīng)、金屬的溶解以及腐蝕產(chǎn)物沉淀。根據(jù)以上問題建立多個(gè)子模型，并將子模型耦合形成LB腐蝕模型。其中，子模型包括電化學(xué)反應(yīng)模型、多組分模型、質(zhì)量傳輸模型、固體邊界遷移模型和結(jié)點(diǎn)鈍化模型。

1.1 電化學(xué)反應(yīng)模型

模擬浸沒于液體腐蝕環(huán)境的金屬（以不銹鋼為例）表面的點(diǎn)蝕現(xiàn)象，其腐蝕過程如圖1所示。

圖1 物理模型示意圖Fig.1 Schematic illustration of corrosion model

為了便于數(shù)值建模，本文將復(fù)雜的電化學(xué)反應(yīng)過程進(jìn)行如下簡化。

陽極：

陰極：

沉淀反應(yīng)：

沉淀氧化：

亞鐵離子氧化：

鐵離子水解：

自催化反應(yīng)：

假設(shè)腐蝕環(huán)境已達(dá)到蝕孔可發(fā)生穩(wěn)態(tài)點(diǎn)蝕時(shí)（腐蝕深度達(dá)到一定值，腐蝕產(chǎn)物溶度達(dá)到其飽和濃度的 60%～80%［8］），如圖1（a）所示。蝕孔內(nèi)金屬基體處于活化溶解狀態(tài)，作為陽極，蝕孔外鈍化膜作為陰極，孔內(nèi)外形成了“大陰極，小陽極”的電化學(xué)腐蝕體系，陰陽極發(fā)生式（1）和式（2）的反應(yīng)，如圖1（b）所示；由于孔外發(fā)生式（2）反應(yīng)，孔口處OH-逐漸增多，孔內(nèi)的金屬離子向孔外擴(kuò)散，發(fā)生式（3）和式（4）的反應(yīng)，在孔口附近產(chǎn)生沉淀；隨著孔內(nèi)金屬不斷腐蝕溶解，金屬離子逐漸增多，繼而發(fā)生式（5）的水解，孔內(nèi)H+增多，呈酸化環(huán)境，進(jìn)一步加快式（6）的腐蝕，如圖1（c）所示。

1.2 多組分模型

使用Shan-Chen模型［2-4］來處理多組分介質(zhì)間的流動(dòng)，其演化方程為

式中：ωα為權(quán)系數(shù)，ω0=4/9，ω1-4=1/9，ω5-8=1/36；宏觀量密度ρ和速度u可由密度分布函數(shù)求得：

在Shan-Chen模型中，為了降低誤差和保證計(jì)算結(jié)果的穩(wěn)定性，使用了分子作用力對速度進(jìn)行修正，分子間的作用力計(jì)算如下：

式中：Fσσ為液體分子間的作用力；Gσσ為作用強(qiáng)度；ψσ為有效密度函數(shù)［9-12］；為權(quán)系數(shù)，其中

式中：Fσw為液體分子與固相分子之間的作用力；wσw為液體與固體壁面的作用強(qiáng)度，可通過接觸角的變化進(jìn)行調(diào)整；s為固體壁面密度［13-15］，與通過分布函數(shù)算出的密度不是同一概念。

將上述計(jì)算得到的分子間作用力代入下式進(jìn)行速度修正：

式中：Fσ為該結(jié)點(diǎn)所受到的合力；u′為式（11）所求的速度；式（14）計(jì)算的u需代入式（9）進(jìn)行平衡分布函數(shù)的計(jì)算。

1.3 質(zhì)量傳輸模型

采用質(zhì)量傳輸模型來解決溶液組分之間的相互傳輸問題，本文涉及到的液體組分為腐蝕溶液組分和腐蝕產(chǎn)物組分。LB質(zhì)量傳輸模型的演化方程為

式中：Ck為第k組分的宏觀量濃度；u由式（11）計(jì)算得到；系數(shù)Jk，α為

其中，Jk，0的取值范圍在0～1之間。

濃度場中的擴(kuò)散系數(shù)Dk為

式中，Jk，0取值為0.9。

1.4 固體邊界遷移模型

由于發(fā)生腐蝕溶解反應(yīng)，被腐蝕的金屬基體表面結(jié)點(diǎn)的體積會(huì)隨之減小，并且當(dāng)新生成的腐蝕產(chǎn)物濃度達(dá)到其飽和濃度時(shí)，在固體表面沉淀結(jié)點(diǎn)的體積會(huì)隨之增大，因此，固體邊界在腐蝕與沉淀的影響下不斷更新。為了能更好地實(shí)時(shí)追蹤新的固體邊界，采用VOP方法予以解決。

根據(jù)式（1）反應(yīng)，金屬基體表面結(jié)點(diǎn)的體積變化［19］如下：

根據(jù)式（2）反應(yīng)，新生成的沉淀結(jié)點(diǎn)在固體表面的體積變化如下：

然后我就發(fā)現(xiàn)，身邊很多女孩跟你一樣，也跟當(dāng)初的我一樣，總是對自己的某一個(gè)缺點(diǎn)耿耿于懷。見人就說：我太胖了，我腿太粗了，我眼睛太小了，我牙齒太丑了……

金屬與腐蝕溶液在邊界處發(fā)生式（1）反應(yīng)的腐蝕，腐蝕邊界處各組分濃度變化［7］為：

式中：CO2(aq)為腐蝕邊界上O2的消耗濃度；為腐蝕邊界上Fe3+的生成濃度；n為反應(yīng)固體表面的外法線方向；為O2擴(kuò)散系數(shù)；為Fe3+擴(kuò)散系數(shù)。

1.5 結(jié)點(diǎn)鈍化模型

模擬過程中的腐蝕邊界，每個(gè)結(jié)點(diǎn)都有一定的概率被鈍化，導(dǎo)致該處的腐蝕速率變慢，使腐蝕坑的形狀呈現(xiàn)不對稱性，這也符合真實(shí)的腐蝕情況。在VOP方法中，腐蝕邊界處結(jié)點(diǎn)發(fā)生體積的變化，如果腐蝕邊界某一結(jié)點(diǎn)發(fā)生鈍化，該結(jié)點(diǎn)的體積減小速率會(huì)變慢。結(jié)點(diǎn)發(fā)生鈍化的概率稱為鈍化概率P，鈍化概率系數(shù)ε為

式中，R為0～1之間的任意隨機(jī)實(shí)數(shù)。因此，需對式（22）進(jìn)行修正：

在VOP方法中引入鈍化概率函數(shù)，可以避免腐蝕界面上的陽極金屬受到的腐蝕作用相同，以使腐蝕坑的形貌隨機(jī)變化且不對稱，更加符合真實(shí)的腐蝕情況。

2 數(shù)值模擬

2.1 物理模型

基于上文構(gòu)建的LB腐蝕模型，對金屬表面點(diǎn)蝕穩(wěn)態(tài)過程進(jìn)行模擬。計(jì)算區(qū)域選取100×160的網(wǎng)格。整個(gè)計(jì)算域內(nèi)，初始的初生蝕孔假設(shè)為一部分圓，其圓心位于（71，80），半徑為10，鈍化膜的厚度為5，如圖1（a）所示。

2.2 邊界條件

對于液體四周的邊界采用了非平衡外推格式進(jìn)行處理，而溶液與金屬之間、溶液與鈍化膜以及溶液與腐蝕產(chǎn)物沉淀之間的流動(dòng)碰撞采用反彈邊界進(jìn)行處理。

2.3 模擬工況

本文模擬金屬表面點(diǎn)蝕穩(wěn)態(tài)過程的形貌的演化過程，參數(shù)如表1所示，表中v為運(yùn)動(dòng)粘度。采用變參數(shù)的數(shù)值模擬得到腐蝕程度的影響因素，各工況參數(shù)如表2所示［3，14］。所有擴(kuò)散系數(shù)均由各濃度場的松弛時(shí)間決定，如式（14）所示。為了滿足運(yùn)行過程中的數(shù)值穩(wěn)定性，選取數(shù)值需大于0.5。

表1 模擬中使用的參數(shù)Table1 Simulation parameters

表2 模擬工況Table 2 Simulated condition

上述工況中所有參數(shù)均為無量綱參數(shù)，均通過表1中的轉(zhuǎn)化關(guān)系得到。由于化學(xué)反應(yīng)速率受溫度、反應(yīng)接觸面積、催化劑等影響，且腐蝕介質(zhì)擴(kuò)散系數(shù)受濃度、溫度等影響，這兩種參數(shù)可選取值范圍較廣［4，15］，因此選取參數(shù)時(shí)，僅在表1中的模擬數(shù)值附近選取。

3 數(shù)值模擬結(jié)果與分析

3.1 腐蝕演化過程分析

通過模擬得到不銹鋼材料表面點(diǎn)蝕穩(wěn)態(tài)過程的形貌變化以及各組分濃度變化情況。模擬工況為k1=8.14×10-5，圖2所示為不銹鋼材料點(diǎn)蝕坑形貌隨時(shí)間變化的情況。從圖中可以看到，腐蝕破損區(qū)域逐漸變大，并且逐漸向金屬基體縱深發(fā)展。當(dāng)蝕孔出現(xiàn)Fe（OH）3沉淀，隨著時(shí)間的推移沉淀增多，逐漸覆蓋了蝕孔（圖2（c））；除了初生的蝕孔發(fā)生腐蝕，在蝕孔底部出現(xiàn)次生蝕孔繼續(xù)向金屬基體深處腐蝕（圖2（d）～圖2（f））。

圖2 腐蝕坑形狀變化Fig.2 Shape evolution of corrosion pit

圖3 腐蝕溶液濃度變化Fig.3 Change in the concentration of corrosion solution

圖3所示為含氧腐蝕溶液的濃度隨時(shí)間變化的情況。隨著時(shí)間的延續(xù)，腐蝕的區(qū)域變大，并向金屬基體縱深發(fā)展，腐蝕溶液向蝕坑內(nèi)滲入，腐蝕溶液逐漸被消耗，在蝕孔處出現(xiàn)沉淀，阻礙了腐蝕溶液向蝕孔內(nèi)擴(kuò)散，導(dǎo)致蝕孔內(nèi)腐蝕溶液的濃度逐漸降低。

圖4所示為蝕孔內(nèi)陽極金屬發(fā)生腐蝕溶解時(shí)，腐蝕產(chǎn)物Fe3+的濃度變化情況。隨著時(shí)間的延續(xù)，腐蝕產(chǎn)物在腐蝕邊界不斷生成，并且向蝕坑外以及液體環(huán)境擴(kuò)散；當(dāng)Fe3+擴(kuò)散至蝕孔與孔外大量的OH-結(jié)合產(chǎn)生Fe（OH）3，F(xiàn)e（OH）3達(dá)到其飽和濃度繼而沉淀。蝕孔產(chǎn)生沉淀，阻礙了Fe3+向孔外擴(kuò)散，孔內(nèi)的Fe3+升高，導(dǎo)致Fe3+進(jìn)一步發(fā)生水解反應(yīng)。

圖4 腐蝕產(chǎn)物Fe3+濃度變化Fig.4 Changes in the concentration of corrosion products Fe3+

圖5所示為蝕孔外鈍化膜作為陰極及其表面產(chǎn)生的OH-濃度變化及OH-擴(kuò)散的情況。隨著時(shí)間的延續(xù)，蝕孔內(nèi)陽極金屬發(fā)生失電子的腐蝕溶解反應(yīng)，電子通過金屬基體轉(zhuǎn)移至電位較高的鈍化膜，O2在鈍化膜表面得到電子，發(fā)生陰極反應(yīng)，在表面產(chǎn)生大量的OH-，并向液體環(huán)境擴(kuò)散；OH-擴(kuò)散至蝕孔處與Fe3+生成Fe（OH）3發(fā)生沉淀反應(yīng)，孔口沉淀的增多阻礙了OH-向蝕孔內(nèi)的擴(kuò)散，導(dǎo)致孔內(nèi)OH-濃度降低。

圖5 OH-濃度變化Fig.5 Changes in the concentration of OH-

圖6所示為H+濃度隨時(shí)間的變化情況。由于蝕孔處的沉淀不斷增多，阻礙了Fe3+向孔外擴(kuò)散以及OH-向孔內(nèi)擴(kuò)散，導(dǎo)致孔內(nèi)大量Fe3+聚集，繼而發(fā)生水解反應(yīng)，產(chǎn)生大量的H+，初生蝕孔底部H+濃度變大，進(jìn)一步腐蝕金屬基體，金屬加速溶解，產(chǎn)生次生孔，推動(dòng)點(diǎn)蝕的穩(wěn)定發(fā)展，蝕孔逐步向深處生長。

圖6 H+濃度變化Fig.6 Changes in the concentration of H+

圖7所示為整個(gè)腐蝕過程中被腐蝕的金屬體積隨時(shí)間的變化關(guān)系。當(dāng)t=0～20 000時(shí)，金屬腐蝕僅發(fā)生初生蝕孔的吸氧腐蝕，初生蝕孔被腐蝕的體積V逐漸變大，蝕孔的深度h逐漸增大；當(dāng)t＞20 000時(shí)，由于孔內(nèi)閉塞及酸化的腐蝕環(huán)境，導(dǎo)致金屬自催化腐蝕反應(yīng)的發(fā)生，腐蝕速率大于先前單獨(dú)發(fā)生的吸氧腐蝕的速率，并且V隨t的增加而增大，h也逐漸增大。

圖7 腐蝕程度隨時(shí)間的變化Fig.7 Variation of the degree of corrosion with respect to time

3.2 影響因素分析

3.2.1 腐蝕化學(xué)反應(yīng)速率

模擬過程中，通過改變式（1）反應(yīng)的化學(xué)反應(yīng)速率，得到某一時(shí)刻（t=50 000）不同腐蝕化學(xué)反應(yīng)速率下點(diǎn)蝕坑的形貌（圖8）。從圖中可以看出，在相同腐蝕時(shí)間下，隨著腐蝕化學(xué)反應(yīng)速率的降低，腐蝕程度減小。

通過對圖8中V和h進(jìn)行統(tǒng)計(jì)測量，得到t=50 000時(shí)V和h與k1的關(guān)系曲線（圖9）。從圖9可以看出，隨著k1增加，V越大，h就越大。由于式（1）中化學(xué)反應(yīng)速率的增加，導(dǎo)致了吸氧腐蝕加快，腐蝕產(chǎn)物增多，蝕孔內(nèi)Fe3+濃度增大，加速了Fe3+水解過程，產(chǎn)生更多的H+，使孔內(nèi)陽極金屬進(jìn)一步加速腐蝕，因此呈現(xiàn)出以上結(jié)果，并且該結(jié)果與真實(shí)情況相符。

圖8 不同腐蝕化學(xué)反應(yīng)速率下的點(diǎn)蝕坑形貌Fig.8 The corrosion pit appearance of different chemical reaction rate

圖9 腐蝕程度隨腐蝕反應(yīng)速率的變化Fig.9 The degree of corrosion varies with corrosion reaction rate

綜合分析圖8與圖9可知，由于k1增加，單位時(shí)間內(nèi)V變大，h增加。

3.2.2 腐蝕溶液擴(kuò)散系數(shù)

通過改變模擬過程中濃度場含氧溶液的擴(kuò)散系數(shù)，研究腐蝕反應(yīng)中介質(zhì)間流動(dòng)與擴(kuò)散對腐蝕的影響。圖10給出了某一時(shí)刻（t=50 000），不同腐蝕溶液擴(kuò)散系數(shù)下腐蝕坑的形貌。從圖中可以看出，相同的腐蝕時(shí)間下，隨著的增大，腐蝕程度加重。

圖10 不同腐蝕溶液擴(kuò)散系數(shù)下的點(diǎn)蝕坑形貌Fig.10 The corrosion pit appearance of different solution diffusion coefficient

通過對圖10中V和h進(jìn)行統(tǒng)計(jì)測量，得到t=50 000時(shí)V和h與DO2的關(guān)系曲線（圖11）。從圖11可以看出，隨著的增大，V越大，h越大。由于的增加，導(dǎo)致腐蝕溶液向蝕孔內(nèi)擴(kuò)散速率加快，金屬腐蝕邊界的濃度相對較高，加快了腐蝕反應(yīng)，因此，呈現(xiàn)出以上模擬結(jié)果。

圖11 腐蝕程度隨O2擴(kuò)散系數(shù)的變化Fig.11 Variation of the degree of corrosion with respect to diffusion coefficient of O2

3.2.3 腐蝕產(chǎn)物擴(kuò)散系數(shù)

在模擬中改變濃度場的腐蝕產(chǎn)物Fe3+的擴(kuò)散系數(shù)，研究腐蝕產(chǎn)物擴(kuò)散系數(shù)對腐蝕程度的影響。圖12所示為某一時(shí)刻（t=50 000）不同腐蝕產(chǎn)物Fe3+擴(kuò)散系數(shù)下腐蝕孔的形貌特征。從圖中可以看出，在相同的腐蝕時(shí)間下，隨著的減小，腐蝕程度加重。

圖12 不同腐蝕產(chǎn)物擴(kuò)散系數(shù)下的腐蝕坑形貌Fig.12 The corrosion pit appearance of different corrosion products diffusion coefficient

通過對圖12中V和h進(jìn)行統(tǒng)計(jì)測量，得到t=50 000時(shí)V和h與的關(guān)系曲線（圖13）。從圖13可以看出，隨著的降低，V越大，h越大。由于的降低，導(dǎo)致蝕孔內(nèi)Fe3+的大量聚集，加速了Fe3+的水解過程，產(chǎn)生了更多的氫離子，使孔內(nèi)的酸度增加，金屬基體的腐蝕溶解速率增加。因此，呈現(xiàn)出以上模擬結(jié)果。

圖13 腐蝕程度隨Fe3+擴(kuò)散系數(shù)的變化Fig.13 Variation of the degree of corrosion varies with respect to diffusion coefficient of Fe3+

4 結(jié) 論

本文基于改進(jìn)后的LB腐蝕模型，將鈍化概率函數(shù)引入VOP方法中，對腐蝕化學(xué)反應(yīng)速率進(jìn)行修正，建立了LB腐蝕模型。使用該模型，模擬了浸沒于液體腐蝕環(huán)境中的金屬表面單坑點(diǎn)蝕穩(wěn)態(tài)過程，通過模擬得到以下結(jié)論：

1）LBM算法可以明確得到溶液中不同組分的流動(dòng)與擴(kuò)散，由此得到其濃度分布，通過流場和濃度場之間的耦合模擬真實(shí)的液體腐蝕環(huán)境；引入鈍化概率函數(shù)使得模擬的腐蝕坑的形貌隨機(jī)變化且不對稱，更加符合真實(shí)的腐蝕坑形貌；陰、陽極共同作用使電化學(xué)腐蝕的模擬更加完善；VOP方法可以控制腐蝕的快慢，而不是簡單的固、液相隨機(jī)轉(zhuǎn)換。

2）腐蝕反應(yīng)速率越大，在相同的腐蝕時(shí)間下，在金屬表面點(diǎn)蝕的穩(wěn)態(tài)過程中，金屬被腐蝕的體積越大，腐蝕的深度就越大；腐蝕含氧溶液擴(kuò)散系數(shù)越大，在相同的腐蝕時(shí)間下，金屬被腐蝕的體積越大，腐蝕的深度就越大；腐蝕產(chǎn)物Fe3+擴(kuò)散系數(shù)越大，在相同的腐蝕時(shí)間內(nèi)，金屬被腐蝕的體積越小，腐蝕的深度就越小。