菱形顆粒沖擊延性材料的運動行為及凹坑形態研究

杜明超， 李增亮， 董祥偉， 孫召成， 范春永， 車家琪

(中國石油大學 (華東)機電工程學院，山東 青島 266580)

微小雜質顆粒(微米至毫米量級)沖擊材料表面導致材料去除的現象廣泛存在于許多工程領域中，例如沖蝕磨損[1]，射流破巖[2]，噴砂清理[3]等。該過程既可以是有益的，也可能是有害的。顆粒沖擊材料表面屬于雙體碰撞問題，其破壞行為取決于沖擊速度、沖擊角度、顆粒形狀及尺寸、材料的類型等因素，可導致材料表面產生微小劃痕(凹坑)、微切削、顆粒破碎、顆粒嵌入等現象。研究者們針對顆粒沖蝕現象進行過許多研究，例如Azimian等[4]建立模擬夾雜顆粒的固液流場環境，結合沖刷實驗研究了沖蝕磨損問題。但是其主要關注一定沖刷時間所造成的沖蝕宏觀效應，并未從顆粒微觀的角度，研究顆粒撞擊微觀過程。

Finnie[5]最早從單個顆粒沖擊入手，研究角型顆粒撞擊導致的材料變形行為，分析了角型顆粒在沖擊過程中對塑性材料的“切削”機制，假定顆粒在碰撞過程中不產生變形，建立了沖蝕理論模型，其適用于顆粒硬度遠高于靶體和低速撞擊的情況。Hutchings等[6-7]改進了Finnie模型，考慮了顆粒在碰撞過程中的自由旋轉，研究了單個方形顆粒沖擊韌性金屬材料表面導致塑性變形特征以及顆粒的反彈行為。同時，Hutchings發現角型顆粒在撞擊過程中的旋轉對顆粒的沖蝕機制有很大的影響。Papini及其研究團隊[8-9]采用氣動發射裝置研究對稱形狀的角型顆粒在射速30～100 m/s的入射速度下的沖蝕機理，并開發了數值計算模型輔助實驗研究，但其研究仍局限在“剛-塑性”假設的框架之內。

近年來，計算流體力學(Computational Fluid Dynamics, CFD)的發展為研究磨料射流、射流破巖以及沖蝕磨損中涉及的固液兩相流問題提供了有效的數值模擬手段。國內外學者開始采用CFD模擬方法對相關過程進行研究。Babets等[10]將水和空氣作為連續相，磨料顆粒作為分散相，采用VOF模型處理連續相，采用拉格朗日(Lagrangian)方法跟蹤顆粒相，得出流場中顆粒軌跡等參數的變化規律。Ahmed等[11]將磨料顆粒、空氣和水都作為連續相，采用Euler方法對其進行處理，借助于商業CFD軟件研究了混和腔內磨料顆粒入射角和入射位置對準直管噴嘴出口處各相速度分布的影響。高激飛等[12]采用雙流體模型，計算了實驗室條件下模擬深水域環境淹沒磨料射流流場，得到了淹沒條件下固液兩相射流的速度與壓力分布規律，并與實驗結果進行了對比。陸朝暉等基于VOF模型和動網格技術建立了截斷式脈沖射流的兩相流瞬態計算模型，研究了射流流體結構的動態演變特征，并與常規圓柱射流進行了對比。基于CFD的數值模擬方法是從固液兩相流場入手，分析固相在流場中的分布，從而得到固相對壁面的沖刷參數，但是再計算沖蝕速率或沖蝕量時，需要借助已有的經驗或半經驗的沖蝕模型。

本文借助ABAQUS商業仿真軟件，建立了菱形顆粒沖擊延性材料的沖蝕磨損耦合數值模型，研究單個菱形顆粒沖蝕磨損過程的運動行為及產生的凹坑輪廓形態，并通過重現參考文獻[9]中的實驗過程驗證了該模型的準確性，隨后針對6組不同沖擊角與方位角的組合研究了沖擊速度對顆粒運動行為及凹坑輪廓的影響，充分了解沖擊速度對沖蝕磨損過程的影響情況，為后續沖蝕機理的研究打下基礎。

1 FEM-SPH模型

1.1 FEM模型

固體力學中的數值模擬方法如有限元法(FEM)可用于顆粒與壁面碰撞過程的模擬。以往的關于微小顆粒沖擊問題的研究大多將顆粒假定為球形，引入變形本構方程描述基底的變形行為，能夠完整重現顆粒沖擊、破壞材料的動態過程，并可實現沖蝕量的定量計算。但FEM并不適合于角型顆粒，這是因為角型顆粒沖擊可造成材料的局部大變形甚至脫落，從而導致網格扭曲等問題，影響數值計算穩定性。此外，FEM迭代時間步長由最小單元尺寸控制，網格極度扭曲會加大計算機CPU耗時，降低工作效率。通過采用自適應網格重構或單元失效等方法可解決上述問題，但自適應網格重構法不能添加失效模型，只適合顆粒前旋沖擊造成的材料堆積情況。其次，計算效率低，每次迭代會重新生成網格，導致整個數值計算過程容易出現錯誤。單元失效法是一種自動刪除扭曲大變形單元的方法，對不同形狀顆粒適用性不同，對尖角菱形顆粒適用性較好，對鈍角方形顆粒適用性差。實驗中顆粒對工件的“犁耕”作用會在工件表面產生材料堆積，仿真模擬中材料堆積處網格單元變形過大被刪除[13]，與物理實際相悖。因此，無網格拉格朗日方法為角型顆粒沖蝕計算提供了另一種解決方案，如SPH方法[14-17]、離散單元法[18](Discrete Element Method, DEM)方法等。

1.2 SPH模型

光滑粒子流體動力學法(Smoothed Particle Hydrodynamics, SPH)是一種無網格(自由網格)的數值計算方法，該方法通過一系列粒子代替網格法中的節點和單元(這些粒子即是計算域中的積分點)，對目標區域進行離散化處理，本質上是一種自適應拉格朗日粒子法。SPH方法中的粒子攜帶所有計算信息，能夠在空間運動，形成了求解描述連續性流體動力學守恒定律的偏微分方程的計算框架。

SPH方法中網格單元被粒子代替，不會出現網格畸變和扭曲纏繞等問題，但邊界區域的處理是SPH方法中一個重要的問題。由于缺少相鄰粒子，標準SPH方程無法求解邊界層粒子的運動方程，為解決該問題， Randles等[19]和Vilad等[20]提出了一種重整化SPH方程，不僅適用于邊界層粒子，也適用于無規則擾動粒子；此外，SPH計算過程耗時長，CPU負荷大，單核運算影響工作效率。

1.3 FEM-SPH模型

由于微小顆粒沖擊的持續時間非常短暫，且時間尺度會隨著顆粒尺寸的降低而降低，在單顆粒實驗中一般使用放大尺寸的顆粒(通常大于1 mm)。因此，有必要引入數值方法來輔助實驗研究更低尺度的顆粒沖擊問題。ABAQUS 6.14是一款功能強大的工程模擬有限元(FE)軟件，在動態顯示求解下，其網格屬性模塊可將網格單元轉化成無網格粒子，用于模擬材料大變形及沖擊載荷的實驗工況。本文采用FEM-SPH耦合方法，采用多處理器并行運算和重整化SPH方程求解，同時解決了FEM和SPH方法的局限性。

顆粒沖擊金屬表面的特征參數定義如圖1所示，vi為初始沖擊速度，αi為沖擊角，θi為初始方位角，θi0為初始旋轉角，vr,αr,θr分別定義為顆粒的反彈速度，反彈角和反彈方位角。A定義為菱形顆粒的菱角度，h為顆粒棱長，w為顆粒厚度，本文中實驗和仿真的數值為：A=60°，h=5.46 mm，w=3.2 mm。

圖1 顆粒特征參數的定義Fig.1 Particle with its characteristic parameters

本文針對Takaffoli和Papini等的實驗對菱形顆粒高速沖蝕過程建立了數值模型，如圖2所示。由圖2可知，在FEM-SPH耦合模型中，只有沖擊后發生材料大變形的區域被離散成無網格粒子，剩余區域以及菱形顆粒仍然離散成有限元網格單元。SPH模塊和FE模塊是通過ABAQUS中tie約束命令耦合在一起。顆粒與金屬表面通過罰定義接觸，為提高仿真結果與實驗結果的吻合度，金屬表面動摩擦因數設為0.1，模型底部完全約束，側面設置關于X-Z平面和Y-Z平面的對稱約束，沖擊顆粒約束為剛體，被沖擊金屬材料為OFHC copper，材料參數，如表1所示。

2 材料模型

2.1 狀態方程(EOS)

固體顆粒沖蝕金屬材料是一個瞬時大變形過程，

圖2 FEM-SPH耦合模型Fig.2 Coupled of FEM-SPH

材料屬性符號數值密度/(kg·m-3)ρ08 960剪切模量/GPaG47.7泊松比υ0.34熔融溫度/KTm1356參考溫度/KTr294比熱/(J·(kgK))CP383J-C初始屈服應力A90 MPaJ-C應變強化系數B292 MPaJ-C硬化指數n0.31J-C應變率強化系數C0.025J-C熱軟化指數m1.09

顆粒沖擊瞬間的動壓遠遠超過材料本身的屈服極限，狀態方程(ESO)通過定義靜壓、密度和比能之間的關系，可完整描述被沖擊材料的動態響應過程[21]。沖擊狀態方程是用于描述固體碰撞的數學方程，需通過動載實驗測出顆粒速度(Up)和沖擊速度(Us)。對于大多數固體和流體而言，在一定壓力范圍內，顆粒速度(Up)和沖擊速度(Us)呈線性關系，但在高速沖擊條件下，顆粒速度(Up)和沖擊速度(Us)的關系表現為非線性，對非金屬材料尤為明顯。ABAQUS中嵌入了Mie-Gruneisen狀態方程，通過輸入OFHC copper材料的Mie-Gruneisen參數可定義被沖擊金屬材料的動態響應過程，OFHC copper材料的Mie-Gruneisen參數，如表2所示。

表2 OFHC copper的Mie-Gruneisen參數Tab.2 M-G parameters of OFHC copper

2.2 Johnson-Cook本構模型

(1)

式中： 無量綱化塑性應變率和溫度定義如公式(2)、式(3)所示。

(2)

T*=T-Tr/Tm-Tr

(3)

式中：ε0是參考等效塑性應變率，Tm和Tr分別為熔融溫度和參考溫度，A，B，C，m，n是材料參數，由實驗獲得，其中參數A，B，n描述了材料屈服應力與應變率和機械硬化之間的關系，參數C和m描述了材料應變率和熱軟化敏感度。本文實驗選用的顆粒材料為W18Cr4V，工件材料為OFHC copper，Johnson-Cook的材料參數，如表1所示。

2.3 失效模型

顆粒后旋沖擊運動對金屬表面產生機加工，出現材料切屑分離和剔除現象，為模擬金屬材料的失效行為，需要引入失效模型。本文采用基于Johnson-Cook塑性模型的失效模型，當方程(4)中破壞參數Di=1[23]時，失效發生，SPH粒子單元i標記失效并從模型表面去除。

D=ΣΔε/εf

(4)

式中： Δε是積分周期內粒子塑性應變的增量，通過計算每個粒子的累積塑性應變增量，得到(ΣΔε)i，εf為材料的失效應變，該失效模型考慮了機械硬化、黏性應變率以及材料熱軟化的影響，其函數方程如公式(5)所示。

(5)

式中：σp為主應力，σe為等效應力，T為材料基體溫度，d1，d2，d3，d4，d5是Johnson-Cook失效模型的5個參數[24]，由實驗測試得出，具體數值見表3。

表3 OFHC copper中Johnson-Cook失效模型參數Tab.3 Johnson-Cook failure model parameters for OFHC copper

3 結果和分析

3.1 模型驗證

為驗證FEM-SPH耦合模型的準確性，首先進行了仿真驗證實驗。表4給出了Takaffoli和Papini等的原始實驗數據，包括初始實驗條件和測得實驗參數。表4中case2、3的顆粒為前旋旋轉(順時針)，反彈角小于90°，顆粒向前反彈；case4、6、7、9、10中顆粒為后旋旋轉(逆時針)，反彈角小于90°，顆粒向前反彈；case1、5中顆粒為后旋旋轉(逆時針)，反彈角大于90°，顆粒向后反彈；case8為特殊工況，顆粒由后旋旋轉轉變為前旋旋轉。保證初始實驗條件不變，進行了與實驗一一對應的仿真模擬，測量出仿真模擬中顆粒的反彈速度(vr)，反彈角度(αr)，凹坑周圍材料堆積高度dmax以及凹坑最大深度hlip，如表5所示。對比表4和表5可知，除case4，其它工況中實驗結果與仿真結果高度吻合。case4中實驗結果出現切屑分離，但仿真模擬并未出現，原因之一是該工況下實際沖擊速度可能大于測量值，隨著仿真模擬中入射速度的增加，切削分離現象出現。

表4 實驗數據Tab.4 Experimental data

表5 仿真模擬預測數據Tab.5 Simulation predicted data

圖3展示了case1、2、8顆粒運動軌跡的實驗結果和仿真結果，實驗結果圖片取自文獻[9]，由高速攝像機抓拍獲得，仿真結果取自本文FEM-SPH耦合模型。對比實驗結果和仿真結果可以看出，該耦合模型能夠很好地捕捉沖蝕磨損過程菱形顆粒的運動行為，3種工況下顆粒運動軌跡的實驗結果與仿真結果吻合度較高，并能夠重現沖蝕磨損過程中顆粒的前旋旋轉和后旋旋轉行為。

圖3 模型驗證Fig.3 Model validation

沖蝕磨損過程顆粒角速度隨時間的變化規律如圖4所示，其中case1、6、10為后旋旋轉，case2為前旋旋轉，case8為后旋向前旋轉變(注：轉變只能從后旋向前旋轉變)。圖5展示了case2和case6中實驗測得的凹坑輪廓與仿真獲得的凹坑輪廓對比，通過對比可以看出仿真模擬結果與實驗測得結果一致，證明了該耦合模型的準確性。

圖4 FEM-SPH模型中顆粒角速度與時間的關系Fig.4 Time history of angular velocities of the particles in FEM-SPH model

3.2 沖擊速度對顆粒運動行為的影響

研究表明，顆粒的沖擊角和方位角是決定顆粒旋向的主要因素，但會受沖擊速度的影響發生輕微改變。沖擊速度對不同旋向顆粒的運動行為影響程度不同。圖6-8展示了前旋旋轉(αi=60°θi=20°，αi=50°θi=20°和αi=40°θi=20°)的顆粒在不同沖擊速度下產生的運動行為以及反彈角速度的變化情況，由圖6～8可知，沖擊速度對前旋顆粒的運動行為和旋轉方向影響較小，顆粒沖擊后的旋轉方向不變，旋轉角速度變化規律一致，符合前旋顆粒沖蝕磨損速率規律[25]。特別需要說明的是，圖6(b)中角速度變化規律表明顆粒接觸到延性材料的表面瞬間有后旋趨勢，隨后轉變為前旋旋轉，這是特殊工況中后旋向前旋轉變的情況，但由于該過程非常短暫，通過實驗很難觀測到初始的后旋過程，且宏觀上顆粒沖擊后表現為向前旋轉，因此劃分為前旋旋轉。

圖5 實驗和仿真預測的case2和case6的凹坑輪廓Fig.5 Experiment and FEM-SPH predicted crater for case 2 and case 6

沖擊速度對“后旋”及臨界沖擊顆粒的運動行為和旋轉方向影響影響較大。圖9～11展示了“后旋”旋轉(αi=60°，θi=50°；αi=60°，θi=50°)和臨界沖擊(αi=40°，θi=50°)的顆粒在不同沖擊速度下產生的運動行為和反彈角速度的變化情況，通過仿真可知，在不同的沖擊速度范圍內，“后旋”顆粒沖擊反彈后的運動行為各不相同，當顆粒入射速度在0～55 m/s時，沖擊后運動行為由后旋向前旋轉變；當顆粒入射速度在55～80 m/s時，顆粒沖擊后出現標準后旋旋轉行為；當入射速度大于80 m/s時，顆粒沖擊后發生前旋旋轉。沖擊速度對該種沖擊角與方位角組合的臨界沖擊影響規律與對“后旋”沖擊的影響規律近似相同。顆粒沖擊后的運動行為和旋轉方向與凹坑形成機制有關，具體原因見下文。

圖6 沖擊速度對前旋旋轉(αi=60°，θi=20°) 顆粒運動行為的影響Fig.6 Effect of impact velocity on kinematicbehavior of particle whenforward rotating(αi=60°，θi=20°)

圖7 沖擊速度對前旋旋轉(αi=50°，θi=20°) 顆粒運動行為的影響Fig.7 Effect of impact velocity on kinematicbehavior of particle whenforward rotating(αi=50°，θi=20°)

圖8 沖擊速度對前旋旋轉(αi=40°，θi=20°) 顆粒運動行為的影響Fig.8 Effect of impact velocity on kinematicbehavior of particle whenforward rotating(αi=40°，θi=20°)

圖9 沖擊速度對后旋旋轉(αi=60°，θi=50°) 顆粒運動行為的影響Fig.9 Effect of impact velocity on kinematicbehavior of particle whenbackward rotating(αi=60°，θi=50°)

圖10 沖擊速度對后旋(αi=50°，θi=50°) 顆粒運動行為的影響Fig.10 Effect of impact velocity on kinematicbehavior of particle whenbackward rotating(αi=50°，θi=50°)

圖11 沖擊速度對臨界沖擊(αi=40°，θi=50°) 顆粒運動行為的影響Fig.11 Effect of impact velocity on kinematicbehavior of particle when critical impact(αi=40°，θi=50°)

3.3 沖擊速度對凹坑輪廓的影響

研究表明[25]，前旋旋轉的顆粒對延性材料表面產生“犁耕”作用，材料堆積在凹坑表面。本文通過數值模擬對三種不同沖擊角與方位角的組合進行研究，獲得了沖擊速度對前旋沖蝕過程產生的凹坑輪廓形態的影響程度，如圖12所示。由圖12可知，同一沖擊角與方位角的組合下，凹坑的尺寸隨著沖擊速度的增長成比例增加，但凹坑的輪廓形態基本不變，這表明，隨著沖擊速度的增加，顆粒尖端對材料產生“犁耕”作用增強，凹坑發生塑性變形的程度更加嚴重，但沖蝕變形機制未發生改變，表明沖擊速度不影響前旋顆粒的沖蝕運動行為。

圖12 前旋沖擊產生的凹坑輪廓Fig.12 Crater profile generated by forward impact

研究表明，“后旋”旋轉的顆粒將導致延性材料“撬起”并切斷剔除，但通過仿真可知，在不同的沖擊速度范圍內，顆粒沖擊獲得的凹坑輪廓形態各不相同。本文通過兩種不同沖擊角與方位角的組合，仿真獲得了沖擊速度對“后旋”沖蝕過程產生凹坑輪廓形態的影響程度，如圖13所示。由圖13可知，當顆粒的入射速度在0～55 m/s時，顆粒初始動能低，不足以將“撬起”的延性材料切斷剔除，未切除的材料對顆粒尖端產生阻礙作用，顆粒運動行為由后旋向前旋轉變；當顆粒的入射速度在55～80 m/s時，延性材料表面發生切屑分離，該過程中顆粒會出現標準后旋旋轉行為，顆粒尖端將“挖掘撬起”的材料切斷剔除，并在材料表面造成一個較淺但較長的凹坑，凹坑輪廓相對比較平滑，深度并不會隨沖擊速度的增加呈明顯增長趨勢；當入射速度在80～100 m/s時，延性材料表面發生切屑分離，顆粒尖端仍然能將“挖掘撬起”的材料切斷剔除，但沖擊后的顆粒發生前旋旋轉，這是因為顆粒初始入射速度較大，法向方向嵌入較深，且顆粒尖端嵌入材料的同時不斷向前堆擠材料，對延性材料表面產生“挖掘”和“微切削”的綜合作用，最終被“挖掘”的材料沿工件表面滑移并切斷，切斷的材料對顆粒尖端產生阻礙作用，導致顆粒發生前旋旋轉；當入射速度大于100 m/s時，顆粒尖端嵌入材料的同時不斷向前堆擠材料，對延性材料表面產生“挖掘”和“微切削”的綜合作用，但初始速度太大，顆粒法向嵌入太深，顆粒尖端的“挖掘力”小于材料抵抗力，被挖掘的材料無法切斷剔除，在“微切削”的作用下產生細長切屑堆積在凹坑前面，導致顆粒發生前旋旋轉。

圖13 “后旋”沖擊產生的凹坑輪廓Fig.13 Crater profile generated by ‘backward’ impact

臨界沖擊的顆粒理論上無前旋或后旋的旋轉趨勢，顆粒沖擊后應按入射路徑返回，但本文研究發現，臨界沖擊的顆粒運動行為和產生的凹坑輪廓受沖擊角與方位角的組合以及沖擊速度的影響較大，仿真證明，對于低沖擊角和大方位角的臨界沖擊組合，不同沖擊速度下顆粒運動行為和凹坑輪廓與“后旋”沖擊產生的規律一致，如圖14所示；對于高沖擊角和小方位角的臨界沖擊組合，不同沖擊速度下顆粒運動行為和凹坑輪廓與前旋沖擊產生的規律一致。由前人的研究成果可知，臨界沖擊下顆粒的動能損失最為嚴重，顆粒旋轉行為被大大抑制，但仍會對接近臨界沖擊的其它沖擊角與方位角的組合存在記憶功能，因此會出現前旋或后旋的沖擊行為。

圖14 臨界沖擊(αi=40°,θi=50°)產生的凹坑輪廓Fig.14 Crater profile generated by critical impact (αi=40°,θi=50°)

顆粒多次沖擊現象是沖蝕磨損過程中一個特殊的現象，顆粒前旋旋轉和后旋旋轉皆會產生多次沖擊，由前期的研究成果可知，在非垂直沖擊和負方位角組合下，前旋顆粒會出現二次乃至三次沖擊，且方位角絕對值越大，顆粒出現的沖擊次數越多；對于顆粒后旋產生切屑分離的工況而言，在非垂直沖擊和正方位角組合下，顆粒至多出現二次沖擊，且方位角越大，二次沖擊效果越明顯，如圖13所示。

4 結 論

(1) 針對菱形顆粒的沖蝕磨損過程建立了基于拉格朗日法的FEM-SPH耦合數值計算模型，通過該模型重現其他研究學者的沖擊實驗，最終實驗測量結果與仿真預測結果吻合度較高，顆粒運動軌跡及反彈參數高度一致，驗證了該耦合模型的準確性。

(2) 沖擊角和方位角是決定顆粒旋轉的關鍵因素，但會受沖擊速度的影響發生變化。沖擊速度對前旋顆粒運動參數的量值影響較大，但對沖擊誘導的顆粒旋轉行為影響較小，不同沖擊速度下導致的凹坑輪廓形狀基本不變，“犁耕”的材料堆積在凹坑表面，符合前旋顆粒沖蝕磨損規律。

(3) 沖擊速度對“后旋”旋轉的顆粒影響較大，當vi在55～100 m/s范圍內時，延性材料表面發生切屑分離，顆粒出現標準后旋旋轉，當vi小于或大于該范圍，只發生材料堆積，均無切屑分離現象發生，低速沖擊產生的材料堆積是由于顆粒初始動能不足無法將“挖掘”的材料切斷剔除，高速沖擊產生的材料堆積是由于顆粒法向嵌入太深，顆粒尖端的“挖掘力”小于材料抵抗力，被挖掘的材料無法切斷剔除。

(4) 臨界沖擊下顆粒的動能損失最為嚴重，顆粒旋轉行為被大大抑制，但仍會對接近臨界沖擊的其它沖擊角與方位角的組合存在記憶功能，因此會出現前旋或后旋的沖擊行為。