不規(guī)則細砂的統(tǒng)計特性及二維/三維數(shù)值模擬

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

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

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

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

目前，國內(nèi)研究沖刷腐蝕的文獻相對較少．劉海笑等[12-13]采用數(shù)值模擬方法分析不同侵蝕公式對彎管侵蝕損傷預測的適用性，研究表明 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]則幾乎全采用宏觀的侵蝕方程進行損傷分布和損傷大小預測，并用形狀因子描述顆粒尖銳程度對侵蝕影響，而未考慮實際顆粒形貌和侵蝕過程．

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

1 試 驗

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

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

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

圖1 SEM檢測結果及相應二值圖Fig.1 Testing result of SEM and the corresponding binary image

顆粒特征參數(shù)統(tǒng)計分布及對應的不同分布函數(shù)擬合優(yōu)度對比結果如圖 2所示．結果表明：對數(shù)正態(tài)(lognormal)和廣義極值(GEV)分布都能較準確地描述砂粒周長p、平面投影面積A、等效直徑D、等效橢圓主軸長Lmaj和短軸長Lmin的分布類型，而縱橫比(AR)則服從GEV分布．該結果與Papini等[9]給出的SiC顆粒檢測結果有略微不同，后者表明粒子的面積服從 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ù)，因此適用于分析復雜分布類型．GEV分布概率密度函數(shù)可表示為

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

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

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

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

圖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ù)相關性Fig.3 Inherent relationship between the distribution parameters of SiO2particles

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

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

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

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

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

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

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

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

圖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ù)預測值(圖 3(b))之間的誤差是否落在周長誤差范圍內(nèi)，若不滿足則返回步驟3．根據(jù)所要模擬的顆粒數(shù)N重復步驟1～5．

圖 7為生成數(shù)值模型的模擬流程．注意：頂點P1、P3、P5和P8在步驟 1確定等效橢圓長、短軸時已完全確定．圖 6(a)中的陰影區(qū)域為頂點控制域(用于象限Ⅰ、Ⅱ、Ⅲ)，象限Ⅳ施加的邊界條件為頂點 3-8、4-8、5-8之間的距離小于主軸長，從而保證顆粒的縱橫比不變．極角，以保證顆粒不出現(xiàn)極其尖銳的頂點(本文取)．試驗檢測結果表明：石英砂顆粒并無特別尖銳的角點，在第 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為試驗結果與數(shù)值模型特征參數(shù)的對比，相應的典型二維數(shù)值模擬結果如圖 9所示．研究表明：數(shù)值模型不僅保證顆粒的縱橫比一致，其顆粒面積、顆粒周長對應的分布與試驗結果吻合良好，表明該數(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檢測結果表明：石英砂顆粒包含片狀、單向平面狀和一般不規(guī)則形狀．假定SiO2顆粒的厚度和 SiC顆粒厚度一樣均服從對數(shù)正態(tài)分布[9]．根據(jù)實際檢測結果，將砂粒分為以下 3種簡化模型：Model Ⅰ：等厚度片狀模型；ModelⅡ：單向變壁厚模型；Model Ⅲ：雙向變壁厚模型．Model Ⅱ和Model Ⅲ的構建是基于Model Ⅰ通過布爾切割運算實現(xiàn)的．根據(jù)第 2節(jié)的數(shù)值模擬算法，采用MCM模擬得到滿足對數(shù)正態(tài)分布的顆粒厚度樣本，通過自定義的 Python程序，采用“自頂向下”的方法通過布爾運算將隨機顆粒模型在ABAQUS中重構，實現(xiàn)厚度不均顆粒的有限元模型構建．三維模型數(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不宜過大，否則會改變顆粒模型的最大厚度，導致砂粒厚度分布不滿足既定分布，而θ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 結 論

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

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

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

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

[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].應用科技，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].應用科技，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] 張 政，程學文，鄭玉貴，等.突擴圓管內(nèi)液固兩相流沖刷腐蝕過程的數(shù)值模擬[J].腐蝕科學與防護技術，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.