王華昆 ，余建星 ，余 楊 ，譚玉娜 ，李修波 ，馮志強

(1.水利工程仿真與安全國家重點實驗室(天津大學(xué))，天津 300072；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240；3.欽州學(xué)院機械與船舶海洋工程學(xué)院，欽州 535011)

管道結(jié)構(gòu)是深海油氣運輸?shù)闹饕绞剑捎谟蜌忾_采中往往混雜大量的酸性氣體(CO2，H2S)和水，形成腐蝕性環(huán)境，此外還夾雜瀝青和固體砂粒，屬于多相流動[1]，在此條件下往往造成嚴重的沖刷腐蝕．沖刷腐蝕過程中，固體顆粒的侵蝕是造成損傷的主要原因．

在含固體顆粒的油氣運輸過程中管道彎管處的沖刷腐蝕尤為嚴重．目前，國外研究結(jié)構(gòu)沖刷腐蝕或侵蝕的主要有 Tulsa大學(xué)的 Vieira等[2-3]、Arabnejad等[4]、Tian 等[5]和 Cheng 等[6]、Takaffoli 等[7-8]及Hadavi等[9-10]．通過將回流式的全管沖刷腐蝕試驗與計算流體動力學(xué) CFD 相結(jié)合，Shirazi等[4]研究了顆粒屬性、流動狀態(tài)、流動介質(zhì)等參數(shù)對彎管沖刷腐蝕的影響．在數(shù)值模擬過程中，假定顆粒為球形顆粒，通過單向或雙向耦合的方式考慮顆粒與流體之間的相互作用．并結(jié)合由試驗得到的經(jīng)驗或半經(jīng)驗侵蝕公式預(yù)測彎管侵蝕損傷及其損傷分布．Tian等[5]和Cheng等[6]采用旋轉(zhuǎn)圓盤電極技術(shù)(RDE)定性研究多相流對沖刷腐蝕的影響．研究表明：油的存在會使鋼表面形成一層油膜，促進了氧的陰極還原并阻礙了鋼的溶解．隨著溶液含砂量和砂粒大小的增大，鋼的腐蝕增強[5]，但文中并未給出砂粒的統(tǒng)計參數(shù)和形貌特征．而 Takaffoli等[7-8]和Hadavi等[9]結(jié)合試驗和數(shù)值模擬重點研究顆粒侵蝕機理和侵蝕過程中靶材的響應(yīng)以及顆粒形貌分析及有限元模擬方法．其中顆粒被簡化為具有均勻厚度的六面體，該數(shù)值模型有明顯的銳角，且不能模擬向內(nèi)凹陷的曲面，與實際不符．Hadavi等[10]根據(jù)試驗結(jié)果考慮了更真實的顆粒形狀對撞擊過程中顆粒破碎的影響．但顆粒輪廓需要通過試驗參數(shù)導(dǎo)入，且模型依然采用等厚度假設(shè).De Pellegrin等[11]開發(fā)了一種模擬具有隨機凸多面體的磨料顆粒方法．這些隨機多面體由平面切割一個立方體得到，平面的方向和位置分布是隨機的．需要指出的是該模型無法模擬含有凹面的多面體，也未將模擬顆粒的統(tǒng)計參數(shù)與實際顆粒統(tǒng)計參數(shù)做對比．

目前，國內(nèi)研究沖刷腐蝕的文獻相對較少．劉海笑等[12-13]采用數(shù)值模擬方法分析不同侵蝕公式對彎管侵蝕損傷預(yù)測的適用性，研究表明 Tabakoff 的侵蝕公式具有較高的精度和較好的適用范圍．盧宇輝等[14-15]則采用CFD模擬研究了通海截止閥的沖刷腐蝕行為．張政等[16]同樣采用數(shù)值模擬方法，考慮了侵蝕過程中的流固雙向耦合作用，分析了突擴圓管在固液兩相流中的沖刷腐蝕．但以上研究均未分析實際的顆粒侵蝕過程，顆粒形狀參數(shù)對侵蝕過程的影響僅采用形狀因子等效考慮．

盡管 Papini等[7-9]對顆粒侵蝕做了較為詳細的微觀統(tǒng)計分析，但用四邊形以及等厚度的板狀三維實體模型模擬 SiC顆粒過于簡化[7-8]，此外管道多相流中攜帶的顆粒多為天然砂粒(SiO2)，其材料屬性與碳化硅(SiC)有較大差異．對于顆粒在侵蝕過程中的損傷分析，Papini等[10]基于試驗測量得到SiC顆粒形貌數(shù)據(jù)，通過拉伸生成厚度均勻的板狀數(shù)值模型，對實際做極大簡化．此外，其他考慮流體作用的沖刷腐蝕研究[2-6,12-16]則幾乎全采用宏觀的侵蝕方程進行損傷分布和損傷大小預(yù)測，并用形狀因子描述顆粒尖銳程度對侵蝕影響，而未考慮實際顆粒形貌和侵蝕過程．

針對已有砂粒數(shù)值模型過于簡化及在只有砂粒分布參數(shù)條件下如何快速建立符合實際的砂粒有限元模型，本文提出新的二維、三維砂粒模型的構(gòu)建方法．基于天然石英砂(SiO2)微粒 SEM 檢測結(jié)果，結(jié)合圖像處理技術(shù)，通過統(tǒng)計分析方法獲取石英砂顆粒特征參數(shù)的分布類型及分布參數(shù)．并通過編程實現(xiàn)簡化顆粒模型的數(shù)值重構(gòu)．通過增加控制點數(shù)量，實現(xiàn)幾何形狀更復(fù)雜的顆粒模擬．此外，新模型考慮了顆粒分布參數(shù)之間的內(nèi)在聯(lián)系．將重構(gòu)模型的特征參數(shù)與試驗數(shù)據(jù)對比驗證模型準確性后，結(jié)合用戶自定義的 Python程序，將不規(guī)則砂粒模型嵌入通用有限元軟件 ABAQUS，實現(xiàn)了不同類型砂粒的有限元模型重構(gòu)，為砂粒沖刷腐蝕數(shù)值模擬研究提供依據(jù)．

1 試 驗

試樣為 6號(40～80目，0.2～0.5mm)天然石英砂顆粒．為建立合理的不規(guī)則石英砂顆粒的有限元模型，首先需獲取顆粒的特征參數(shù)統(tǒng)計分布．將石英砂顆粒表面噴金處理后，采用 SEM 進行檢測，檢測結(jié)果如圖1(a)所示．

為便于對石英砂顆粒進行形貌分析，首先將SEM 掃描圖片轉(zhuǎn)化為二值圖，如圖 1(b)所示，并用Photoshop提取顆粒輪廓(僅提取輪廓清晰的圖像).根據(jù) SEM 掃描圖對應(yīng)的標尺，將所有樣本處理成1mm×1mm的黑底圖片，每個圖片包含一個顆粒，共 418個樣本．通過 Matlab編程逐個提取顆粒形貌信息并存儲，最后匯總進行統(tǒng)計分析．其中面積參數(shù)采用“regionprops”函數(shù)提取，等效直徑(D)和顆?？v橫比(aspect ratio，AR)定義為

式中：A為石英砂顆粒平面投影面積；Lmaj為等效橢圓主軸長；Lmin為等效橢圓短軸長．等效橢圓定義為與顆粒平面投影形狀具有相同二階矩的橢圓．

圖1 SEM檢測結(jié)果及相應(yīng)二值圖Fig.1 Testing result of SEM and the corresponding binary image

顆粒特征參數(shù)統(tǒng)計分布及對應(yīng)的不同分布函數(shù)擬合優(yōu)度對比結(jié)果如圖 2所示．結(jié)果表明：對數(shù)正態(tài)(lognormal)和廣義極值(GEV)分布都能較準確地描述砂粒周長p、平面投影面積A、等效直徑D、等效橢圓主軸長Lmaj和短軸長Lmin的分布類型，而縱橫比(AR)則服從GEV分布．該結(jié)果與Papini等[9]給出的SiC顆粒檢測結(jié)果有略微不同，后者表明粒子的面積服從 Gaussian分布，而D服從 Weibull分布[7-9]．實際上 Weibull分布是 GEV分布的特殊情況，因此D也服從GEV分布．采用Gaussian分布描述顆粒面積的主要不足是：在蒙特卡洛(MCM)模擬過程中，模擬參數(shù)可能會出現(xiàn)負值，這與實際不符．因此，在模擬過程中往往需要對模擬值符號進行判斷與限界，而采用對數(shù)正態(tài)分布可避免該問題．對數(shù)正態(tài)分布的概率密度函數(shù)為

其分布均值m和方差v可通過分布參數(shù)μ和σ得到.

由于GEV分布有3個分布參數(shù)，因此適用于分析復(fù)雜分布類型．GEV分布概率密度函數(shù)可表示為

式中：μ為位置參數(shù)；σ為尺度參數(shù)；k為形狀參數(shù)．

對于，當k＞0，對應(yīng)Ⅱ型廣義極值分布；當k＜0，對應(yīng)Ⅲ型廣義極值分布；當k＝0，對應(yīng)Ⅰ型分布，即極值分布，此時對應(yīng)概率密度函數(shù)為

由于GEV分布具有廣泛適用性，因此采用GEV分布對各分布參數(shù)的試驗統(tǒng)計數(shù)據(jù)進行擬合．對應(yīng)的分布參數(shù)如表1所示．

實際上顆粒分布參數(shù)之間并非相互獨立，面積A和周長p都依賴于顆粒的長、短軸．根據(jù)試驗檢測結(jié)果，通過非線性回歸分析得到A和p與顆粒的長、短軸之間的關(guān)系如式(6)所示．解析函數(shù)與試驗對比結(jié)果及對應(yīng)的擬合優(yōu)度參數(shù)見圖3．圖4給出了解析公式預(yù)測值與試驗結(jié)果對比的誤差分布及采用GEV分布擬合得到的回歸參數(shù)．結(jié)果表明：該經(jīng)驗公式預(yù)測值與試驗結(jié)果吻合良好，大部分誤差均在±5%以內(nèi)，且GEV分布可準確擬合誤差分布．圖5給出了考慮擬合誤差和不考慮擬合誤差的模擬結(jié)果對比．結(jié)果表明二者預(yù)測值差異很小．

圖2 SiO2顆粒特征參數(shù)統(tǒng)計分布Fig.2 Statistical distribution of particle characteristic parameters of SiO2particles

表1 石英砂顆粒特征分布參數(shù)Tab.1 Characteristic distribution parameters of the silica sand particles

圖3 SiO2顆粒特征參數(shù)相關(guān)性Fig.3 Inherent relationship between the distribution parameters of SiO2particles

在 Takaffoli等[7]的研究中，并未考慮顆粒面積、周長與顆粒長短軸之間的內(nèi)在聯(lián)系，也未考慮顆?？v橫比應(yīng)滿足的關(guān)系．僅通過球度參數(shù)(S＝4πA/p2)考慮面積、面積與周長之間的關(guān)系，這種約束是偏弱的．實際上縱橫比是顆粒的主要特征，該參數(shù)決定了顆粒是條狀還是近似為方形．針對以上不足，本文模型通過考慮分布參數(shù)之間的內(nèi)在聯(lián)系，使得顆粒形狀更加合理．并通過增加控制點數(shù)量實現(xiàn)更復(fù)雜幾何形狀的模擬．數(shù)值模擬方法見第2節(jié)．

圖4 回歸分析誤差分布Fig.4 Error distribution of regression analyis

圖5 考慮和不考慮擬合誤差的模擬結(jié)果對比Fig.5 Comparison of the simulation results with and without the consideration of regression error

2 二維數(shù)值模擬與驗證

根據(jù)第1節(jié)試驗統(tǒng)計分析結(jié)果，采用Matlab編程實現(xiàn)平面顆粒模型二維重構(gòu)．采用8頂點模型模擬顆粒平面投影形狀，并根據(jù)顆粒厚度分布隨機數(shù)控制三維簡化模型(見第 3節(jié))．二維數(shù)值模型模擬過程如下.

步驟 1 根據(jù)表 1主軸長和縱橫比滿足的分布類型及對應(yīng)的分布參數(shù)，采用 MCM 隨機生成樣本Lmaj、AR，從而確定短軸長Lmin＝Lmaj/AR；

步驟 2 結(jié)合式(6)的回歸分析結(jié)果，根據(jù)主軸、短軸長確定顆粒面積(A0)，考慮回歸分析的誤差分布(圖 4(a))，采用 MCM 得到面積預(yù)測誤差樣本ErrorA，從而確定實際用于模擬的顆粒面積Ac＝A0(1＋ErrorA)；

步驟 3 假設(shè)極徑Ri和極角θi均滿足均勻分布，采用MCM隨機生成樣本Ri和θi(圖6(a))，從而確定象限Ⅰ～Ⅳ4個頂點的X、Y坐標(Pi表示頂點P編號，見圖 6(b))，各參數(shù)計算方法及對應(yīng)參數(shù)范圍如式(7)所示．4個頂點坐標信息均在象限Ⅰ中生成，然后變換至其他象限．θmin表示控制點P與原點O構(gòu)成的直線與橫軸之間夾角的最小容許值．

圖6 數(shù)值模型示意Fig.6 Schematic of the numerical model

步驟 4 根據(jù)顆粒面積Ac、8個頂點(P1～P8)的X坐標以及前7個頂點(P1～P7)的Y坐標，重置第8個頂點的Y坐標．

步驟5 計算顆粒頂點3-8間距(L38)、4-8(L48)、5-8(L58)及顆粒周長，判斷實際周長與回歸函數(shù)預(yù)測值(圖 3(b))之間的誤差是否落在周長誤差范圍內(nèi)，若不滿足則返回步驟3．根據(jù)所要模擬的顆粒數(shù)N重復(fù)步驟1～5．

圖 7為生成數(shù)值模型的模擬流程．注意：頂點P1、P3、P5和P8在步驟 1確定等效橢圓長、短軸時已完全確定．圖 6(a)中的陰影區(qū)域為頂點控制域(用于象限Ⅰ、Ⅱ、Ⅲ)，象限Ⅳ施加的邊界條件為頂點 3-8、4-8、5-8之間的距離小于主軸長，從而保證顆粒的縱橫比不變．極角，以保證顆粒不出現(xiàn)極其尖銳的頂點(本文取)．試驗檢測結(jié)果表明：石英砂顆粒并無特別尖銳的角點，在第 3節(jié)有限元模型中將通過增加圓角進一步改進．

ABAQUS 的內(nèi)核語言是 Python語言，由于Python中沒有可生成滿足 GEV分布的內(nèi)置隨機函數(shù)，因此，通過求累積分布函數(shù)F(x)的逆函數(shù)得到縱橫比樣本x．此外，也可直接將Matlab生成的幾何參數(shù)直接寫入.py文件，進行參數(shù)化建模．給定頻率y的GEV分布逆函數(shù)可表示為

圖7 二維模型數(shù)值模擬流程Fig.7 Flow chart of the numerical simulation of two-dimensional particle

圖8為試驗結(jié)果與數(shù)值模型特征參數(shù)的對比，相應(yīng)的典型二維數(shù)值模擬結(jié)果如圖 9所示．研究表明：數(shù)值模型不僅保證顆粒的縱橫比一致，其顆粒面積、顆粒周長對應(yīng)的分布與試驗結(jié)果吻合良好，表明該數(shù)值算法是可靠的．第 3節(jié)將基于該二維平面模型建立三維有限元模型．

圖8 2D模型驗證Fig.8 Validation of two-dimensional model

圖9 數(shù)值模擬二維粒子幾何形狀Fig.9 Simulation results of the two-dimensional particles

3 三維模型有限元模擬

二維模型常常作為三維模型的簡化，第2節(jié)通過數(shù)值模擬已得到合理的不規(guī)則砂粒平面投影二維模型．顆粒形狀是影響侵蝕過程的主要因素之一，不同的顆粒形狀引起的侵蝕機理也不同：對于尖銳粒子來說，主要觀察到韌性侵蝕，而在球狀顆粒中觀察到脆性侵蝕．此外，有角顆粒會嵌入基體材料，而球形粒子一般不會[17]．由于砂粒形狀會影響侵蝕機理從而影響侵蝕過程，因此，建立砂粒的三維隨機模型對砂粒侵蝕機理研究具有重要的意義．

圖1的顆粒形貌SEM檢測結(jié)果表明：石英砂顆粒包含片狀、單向平面狀和一般不規(guī)則形狀．假定SiO2顆粒的厚度和 SiC顆粒厚度一樣均服從對數(shù)正態(tài)分布[9]．根據(jù)實際檢測結(jié)果，將砂粒分為以下 3種簡化模型：Model Ⅰ：等厚度片狀模型；ModelⅡ：單向變壁厚模型；Model Ⅲ：雙向變壁厚模型．Model Ⅱ和Model Ⅲ的構(gòu)建是基于Model Ⅰ通過布爾切割運算實現(xiàn)的．根據(jù)第 2節(jié)的數(shù)值模擬算法，采用MCM模擬得到滿足對數(shù)正態(tài)分布的顆粒厚度樣本，通過自定義的 Python程序，采用“自頂向下”的方法通過布爾運算將隨機顆粒模型在ABAQUS中重構(gòu)，實現(xiàn)厚度不均顆粒的有限元模型構(gòu)建．三維模型數(shù)值模擬流程如圖10所示．其中：切割部件與某一邊平行，夾角θj(θmin≤θj≤θmax)定義為切割面與豎軸之間的夾角，如圖 11所示，圖中陰影區(qū)為切割區(qū)域．

圖10 三維模型數(shù)值模擬流程Fig.10 Flow chart of the numerical simulation of three-dimensional model

根據(jù)類型索引號(圖9中i的值)按圖10流程圖生成不同類型的三維微粒模型．其中θmax不宜過大，否則會改變顆粒模型的最大厚度，導(dǎo)致砂粒厚度分布不滿足既定分布，而θmin可不做限制，θmin＝0表示不切割．切割深度h1和h2是隨機確定的，計算方式為

式中：h為顆粒厚度且滿足對數(shù)正態(tài)分布分布；Llow和Lup分別表示隨機數(shù)下限和上限，一般取0和1．

圖11 布爾切割運算示意Fig.11 Schematic of Boolean cutting operation

圖 12(a)～(c)分別給出了 3種簡化模型的幾何模型和有限元模型．從左至右依次為一般模型、圓角模型和有限元模型．圖 12(d)中給出的為全部施加圓角后的模型，在模擬過程中可根據(jù)模型幾何邊索引號(隨機確定)對任意幾何邊進行不同圓角半徑過渡，使數(shù)值模型更貼近實際．圖 12(b)為單向變壁厚模型，圖 12(c)為雙向變壁厚模型，圖 12(d)為典型的不同類型顆粒有限元模型．從圖 12中可以看出該方法可得到任意幾何形狀的顆粒，后續(xù)將基于該顆粒模型展開深入的顆粒侵蝕機理分析．

圖12 細砂顆粒三維數(shù)值模型Fig.12 Three-dimensional numerical model of fine sand particle

4 結(jié) 論

本文基于試驗檢測結(jié)果得到天然石英砂微粒(40～80目)特征參數(shù)的分布類型和對應(yīng)的分布參數(shù)．根據(jù)試驗結(jié)果提出一種新的二維、三維不規(guī)則顆粒模型模擬方法，得到如下結(jié)論．

(1) 砂粒周長p、平面投影面積A、等效直徑D、等效橢圓主軸長Lmaj和短軸長Lmin服從 Lognormal或GEV分布，而縱橫比AR則服從GEV分布．GEV分布具有最廣泛的適用性，適用于描述復(fù)雜分布．

(2) 新數(shù)值方法通過考慮顆粒長短軸之間的比例關(guān)系以及面積、周長與長短軸之間的內(nèi)在聯(lián)系，使得數(shù)值模型更貼近實際．

(3) 基于二維模型建立的三維模型既保證平面特征參數(shù)的統(tǒng)計特性，又考慮了顆粒厚度分布．通過不同的簡化模型，編程實現(xiàn)了不同類型石英砂微粒三維數(shù)值重構(gòu)，通過圓角功能避免了不切實際的銳角，為管道的顆粒侵蝕損傷機理研究提供必要的信息．

[1] Ige O O，Umoru L E.Effects of shear stress on the erosion-corrosion behaviour of X-65 carbon steel：A combined mass-loss and profilometry study[J].Tribology International，2016，94(2)：155-164.

[2] Vieira R E，Parsi Mazdak，Zahedi Peyman，et al.Sand erosion measurements under multiphase annular flow conditions in a horizontal-horizontal elbow[J].Powder Technology，2017，320：625-636.

[3] Vieira R E，Mansouri A，Mclaury B S，et al.Experimental and computational study of erosion in elbows due to sand particles in air flow[J].Powder Technology，2016，288：339-353.

[4] Arabnejad H，Shirazi S A，Mclaury B S，et al.The effect of erodent particle hardness on the erosion of stainless steel[J].Wear，2015，332/333：1098-1103.

[5] Tian B R，Cheng Y F.Electrochemical corrosion behavior of X-65 steel in the simulated oil sand slurry.Part I：Effects of hydrodynamic condition[J].Corrosion Science，2008，50(3)：773-779.

[6] Tang X，Xu L Y，Cheng Y F.Electrochemical corrosion behavior of X-65 steel in the simulated oil-sand slurry.PartⅡ：Synergism of erosion and corrosion[J].Corrosion Science，2008，50(5)：1469-1474.

[7] Takaffoli M，Papini M.Numerical simulation of solid particle impacts on Al6061-T6，Part I：Threedimensional representation of angular particles[J].Wear，2012，292/293(29)：100-110.

[8] Takaffoli M，Papini M.Numerical simulation of solid particle impacts on Al6061-T6，Part Ⅱ：Materials removal mechanisms for impact of multiple angular particles[J].Wear，2012，296(1/2)：648-655.

[9] Hadavi V，Michaelsen B，Papini M.Measurements and modeling of instantaneous particle orientation within abrasive air jets and implications for particle embedding[J].Wear，2015(336/337)：9-20.

[10] Hadavi V，Moreno C E，Papini M.Numerical and experimental analysis of particle fracture during solid particle erosion，Part I：Modeling and experimental verification[J].Wear，2016(356/357)：135-145.

[11] De Pellegrin D V，Stachowiak G W.Simulation of threedimensional abrasive particles[J].Wear，2005，258：208-216.

[12] 王思邈，劉海笑，張 日，等.海底管道沙粒侵蝕的數(shù)值模擬及侵蝕公式評價[J].海洋工程，2014，32(1)：49-59.

Wang Simiao，Liu Haixiao，Zhang Ri，et al.Numerical simulations of sand erosion in pipelines and evaluations of solid particle erosion equations[J].The Ocean Engineering，2014，32(1)：49-59(in Chinese).

[13] 張 日，劉海笑.流動保障中管道的顆粒侵蝕分析[J].海洋工程，2012，30(4)：10-20.

Zhang Ri，Liu Haixiao.Solid particle erosion analysis of pipelines in flow assurance[J].The Ocean Engineering，2012，30(4)：10-20(in Chinese).

[14] 盧宇輝，盧銳新，程千駒，等.影響通海截止閥沖刷腐蝕的重要參數(shù)數(shù)值模擬分析[J].應(yīng)用科技，2016，43(4)：70-75，80.

Lu Yuhui，Lu Ruixin，Cheng Qianju，et al.Numerical simulation and analysis of key parameters influencing the scouring corrosion of the deep diving cutoff valve[J].Applied Science and Technology.2016，43(4)：70-75，80(in Chinese).

[15] 盧宇輝，盧銳新，劉少剛.通海截止閥沖刷腐蝕的數(shù)值模擬及優(yōu)化[J].應(yīng)用科技，2016，43(1)：61-66.

Lu Yuhui，Lu Ruixin，Liu Shaogang.Numerical simulation and optimization for erosion corrosion of the cutoff valve[J].Applied Science and Technology，2016，43(1)：61-66(in Chinese).

[16] 張 政，程學(xué)文，鄭玉貴，等.突擴圓管內(nèi)液固兩相流沖刷腐蝕過程的數(shù)值模擬[J].腐蝕科學(xué)與防護技術(shù)，2001(2)：89-95.

Zhang Zheng，Cheng Xuewen，Zheng Yugui，et al.Numerical simulation of errosion-corrosion in liquidsolid two-phase flow[J].Corrosion Science and Protection Technology，2001(2)：89-95(in Chinese).

[17] Roy M，Tirupataiah Y，Sundararajan G.Effect of particle shape on the erosion of Cu and its alloys[J].Materials Science and Engineering：A，1993，165(1)：51-63.