磁致伸縮儀超聲空化流的三維非定常數值模擬

夏冬生 孫昌國 于 彥 張會臣

大連海事大學，大連，116026

?

磁致伸縮儀超聲空化流的三維非定常數值模擬

夏冬生 孫昌國 于 彥 張會臣

大連海事大學，大連，116026

采用Singhal完全空化模型和SSTk-ω湍流模型結合動網格技術對磁致伸縮儀超聲空化流進行三維非定常數值模擬。計算結果表明，由于變幅桿高頻振動，在靠近試樣表面附近局部流場的壓力和空泡體積組分變化具有周期性，壓力波動的最低值可達到汽化壓力，該局部流場可發生空化。由于空化，試樣表面壓力波動具有脈沖特征。壓力和空泡體積組分在試樣表面近似呈環形分布。在同一環形區域內，壓力和空泡體積組分存在無規律斷續脈動。試樣表面中心區域空泡經歷兩次振蕩后潰滅，產生強烈脈沖壓力，最大脈沖壓力可達約14 MPa。脈沖壓力在試樣表面按間隔環形區域分布，且隨試樣振動在相鄰環形區域上交替出現。在磁致伸縮儀超聲空化流場中，試件表面可近似多個聲波發生源，各聲波傳播時相互疊加和干擾。在聲波傳播的過程中壓力衰減很快，只是在距試樣表面約20 mm內，壓力有明顯波動。

磁致伸縮儀；計算流體動力學；空化；脈沖壓力

0 引言

空蝕是流體機械的一種主要破壞形式，嚴重影響流體機械的工作可靠性、工作效率和使用壽命。國內外學者對空蝕進行了廣泛研究，其研究領域包括空蝕機理[1]、空化流場特征[2]、材料抗空蝕性能[3]及空化和空蝕影響因素[4]等。目前，研究者普遍接受的空蝕機制是微射流理論[5-6]和沖擊波理論[7]。當空泡在固體表面附近發生潰滅時，產生高速微射流或(和)沖擊波，其作用在固體表面形成顯著的沖擊壓力。表面材料由于受到沖擊壓力的頻繁作用而發生疲勞破壞。空化流流場特性和材料的抗空蝕能力決定了空蝕破壞是否發生及破壞的程度。

磁致伸縮儀是研究空蝕機理和測定材料抗空蝕性能的重要儀器之一。磁致伸縮儀的基本原理是利用具有趨磁性的傳感器或壓電傳感器在交變電流作用下能夠伸長或變短的特性，實現換能器端部在液體中高頻振動，產生振蕩型無主流空化。為了深入理解磁致伸縮儀空蝕實驗中材料的空蝕行為及空蝕特征，研究振蕩型空化流空化過程和空化流流場特征是十分必要的。Ahmed[4]利用磁致伸縮儀研究了溫度對沖擊壓力、空化過程和金屬鋁空蝕行為的影響，發現空泡潰滅引起的沖擊壓力分布在試樣表面圓形區域內，且與溫度和徑向距離密切相關：溫度較低時，中心區域沖擊壓力最大；隨溫度升高，中心區域的沖擊壓力減小，易形成較多數量的空泡。Bai等[8]利用磁致伸縮儀和高速攝影技術對超聲空化流流場中硬壁凹坑的駐氣特性、微泡散逸特性以及空泡對空蝕坑發展過程的影響進行了研究。Burdin等[9]利用磁致伸縮儀和激光技術研究了超聲空化的空泡體積含量、空泡大小分布以及超聲功率的影響。由于空化的瞬時性、隨機性和復雜性以及實驗條件限制，磁致伸縮儀超聲空化過程和流場特征的研究較缺乏。本文基于CFD方法，利用Fluent軟件模擬磁致伸縮儀超聲空化流流場，分析超聲空化流場的壓力、空泡分布等流場細節和空化發展過程，并根據空化流的流場特性來探討磁致伸縮儀空蝕實驗中試樣的空蝕基本特征。

1 數學模型

1.1 空化流的控制方程

空化流一般由液體和含有蒸氣、不可凝結氣體的空泡群組成。本文采用基于均質平衡流理論的Singhal完全空化模型模擬超聲空化流。該空化模型考慮了不可凝結氣體的影響，將空化流視為由液體、蒸氣和不可凝結氣體三組元組成的氣-液兩相混合的單流體(即兩相混合的均勻介質)。三組元共享同一物理場，如壓力、速度等。假設空化過程為等溫過程，忽略重力效應，則混合單流體的Favre平均的連續方程和動量方程分別為

(1)

(2)

其中，xi、xj、xk(i,j,k=1, 2, 3)為笛卡兒坐標，下標1、2和3分別代表x軸、y軸和z軸方向，ui是沿i方向的速度分量，p是壓力，μ、μt分別為混合流體的動力黏度和湍流黏度，δij是克羅內克(Kronecker)數。

混合流體的密度ρm由下式計算：

(3)

其中，ρ是密度，f是質量分數，下標l、v和g分別代表液體、蒸氣和不可凝結氣體。

不可凝結氣體一般選擇空氣并假設滿足理想氣體狀態方程，其密度ρg由公式ρg=Wp/(RT)確定，其中，T為溫度，通用氣體常數R為8314.3 J/(kg·mol)·K，W為不可凝結氣體的分子量。

空化過程中液-氣相間的質量傳輸由蒸氣質量分數fv的輸運方程控制：

(4)

式中，Re、Rc分別為蒸發階段(蒸氣泡的產生和膨脹)和凝結階段(蒸氣泡的壓縮和破裂)的源項，分別定義為蒸氣的蒸發率和凝結率。

對于Singhal完全空化模型，Re和Rc從Rayleigh-Plesset方程推導得到，考慮了湍動能在空泡運動過程中的影響和不可凝結氣體作為氣核的影響等，其表達式如下[10]：

(5)

(6)

其中，σ為液體的表面張力系數；k為流場的當地湍動能；經驗常數Ce=0.02，Cc=0.01。Singhal完全空化模型考慮了湍流壓力脈動的影響，修正相變壓力閥值pv從飽和蒸氣壓psat提高到pv=psat+pturb/2，其中，pturb為湍流引起的壓力脈動，pturb=0.39ρmk。

1.2 剪切應力輸運(SST)k-ω湍流模型

采用SSTk-ω湍流模型對以上控制方程封閉。SSTk-ω湍流模型是由k-ω模型和k-ε模型各乘以一混合函數相加而成的，其在近壁區等價于低雷諾數標準k-ω模型，在遠離壁面區域等價于高雷諾數k-ε模型，而在混合區域由混合函數混合使用這兩種模型[11]。SSTk-ω模型除了兼具以上兩種模型在近壁區和遠場計算的優點外，與k-ω模型相比，ω輸運方程中還增加了交叉擴散項，同時在湍流黏度定義中考慮了湍流剪切應力的輸運過程，從而使該模型適用范圍更廣。

SSTk-ω模型中湍動能k和比耗散率ω的輸運方程分別為

(7)

(8)

(9)

其中，S為平均應變率的張量模量，α、β*、β均為模型系數。湍流黏度μt按下式計算：

(10)

其中，對于高雷諾數，α*=1；對于低雷諾數，α*被修正，模型常數a1=0.31。

k和ω的湍流普朗特數Prk和Prω分別為

(11)

(12)

2 計算模型、網格劃分和計算方法

2.1 計算模型和邊界條件

圖1所示為磁致伸縮儀空蝕實驗裝置，測試采用試樣固定在變幅桿端面隨變幅桿振動的方式。磁致伸縮儀的典型工作參數如下：變幅桿振動頻率為20 kHz，峰-峰振幅為50～75 μm，功率為1 kW。忽略機電控制單元和冷卻恒溫裝置，建立磁致伸縮儀空化流的計算模型，其模型尺寸如圖2所示。變幅桿端部為直徑20 mm的圓平面，浸入水深為20 mm。為了使聲場近似為自由場，減少或避免聲波的反射影響，將容器的直徑和水深均設置為變幅桿直徑的5倍。

圖1 磁致伸縮儀空蝕實驗裝置示意圖

圖2 三維計算模型

定義自由液面為壓力進口，壓力設定為101 kPa。所有壁面均采用無滑移邊界條件。

利用動網格技術模擬變幅桿端面的振動，變幅桿端面的運動規律由用戶自定義函數(user’s defined function，UDF)來描述。由于變幅桿端面運動為剛體平動，故調用DEFINE_CG_MOTION宏函數來編寫UDF程序文件，需要指定變幅桿端面振動的速度。

根據變幅桿端面振動的運動規律，同時為有利于變幅桿振動模擬的實現，變幅桿端面的往復運動采用簡諧運動近似表示，其位移方程為

z=Asin(2πft)

(13)

其中，z為振動位移；A為振幅；f為振動頻率；t為時間。根據磁致伸縮儀的工作參數，計算時選取振動頻率f為20 kHz，振幅A為30 μm。

將式(13)兩邊對時間求導，得到變幅桿端面的運動速度方程：

v=2πfAcos(2πft)

(14)

將編寫好的UDF程序源文件調入Fluent中進行編譯和連接，實現變幅桿端面按UDF文件所定義的方式運動。變幅桿端面的速度和加速度隨時間變化曲線如圖3所示。

2.2 網格劃分

為提高數值模擬準確性和計算精度，對計算域采用六面體結構化網格劃分，并對變幅桿端面和側面附近網格進行加密處理。由于變幅桿振幅為30 μm，同時SSTk-ω模型要求壁面邊界層網格的量綱一壁面距離y+在較低范圍，則變幅桿端面邊界層網格的首層高度設定為20 μm。但需要注意的是，由于變幅桿端面運動，其邊界層網格發生畸變，將使變幅桿端面y+值發生顯著變化，為解決這一問題，對所有網格類型采用彈簧近似光順法來更新網格，并設置彈簧彈性系數為0，實現動邊界附近的網格保持其原始形狀和密度隨動邊界一起移動，網格變形在距動邊界一定區域外進行。計算結果表明，變幅桿端面處第一層網格節點的y+值在0.7～9.0之間，滿足SSTk-ω模型對y+的要求。流場計算域的網格劃分如圖4所示，整個計算域共得到約330萬個網格。

(a)振動速度變化曲線

(b)振動加速度變化曲線圖3 變幅桿端面的振動速度和加速度隨時間變化曲線

(a)縱截面(b) 過試樣表面的橫截面圖4 流場計算域網格

2.3 計算方法

采用有限體積法對控制方程組進行離散，對流項的離散采用二階迎風格式，壓力離散方式采用PRESTO!，擴散項的離散采用具有二階精度的中心差分格式，時間離散采用一階全隱式格式。為了提高非定常空化流計算的收斂速度和計算穩定性，先不考慮變幅桿振動進行定常計算，其收斂結果作為非定常計算的初始流場。定常和非定常計算分別采用SIMPLEC算法和PISO算法實現速度和壓力之間的耦合求解。考慮到超聲空化流場特性是瞬時變化的，為提高計算準確性和捕捉流場變化細節，本文將時間步長設定為變幅桿振動周期T的1%，即0.5 μs。本文對變幅桿振動20個振動周期內空化流場進行數值計算。

計算中空化流的物性參數取溫度為30 ℃時的值，其中水和水蒸汽的密度分別為995.6 kg/m3和30.36 g/m3，水和水蒸汽的動力黏度分別為8×10-4kg/(m·s)和1.34×10-5kg/(m·s)，水的飽和蒸汽壓psat為4247 Pa，水的表面張力系數σ為0.0717 N/m。假定水中不可凝結氣體質量分數fg為1.5×10-5(默認值)。

3 計算結果與分析

在試件表面分別選擇中心、邊緣和徑向距離5 mm處三個特征點，監視其壓力和氣體體積組分(vapor volumn fraction, VVF)隨時間波動情況，其變化曲線分別如圖5和圖6所示。從圖5和圖6可看出，隨變幅桿振動，試件表面各處壓力和氣體體積含量的波動具有周期性，而且壓力波動具有顯著的脈沖特征。除試件邊緣外，其他兩處壓力可減小到水的汽化壓力附近，表明試件表面附近流場可發生空化。試件表面附近流場空化表現為瞬態空化的特征。在試件中心處隨壓力的波動，空泡經歷一次膨脹、收縮的振蕩過程后，在超聲的負壓力區第二次膨脹，空泡體積進一步增加，VVF值達到50%，直至在隨后出現的正壓力區潰滅。空泡潰滅時，該處產生脈沖壓力，峰值可達14 MPa。試樣中心的脈沖壓力每兩個振動周期出現一次，這與該處的空泡潰滅頻次相一致。在徑向距離為5 mm處，一個振動周期內大部分空泡只經歷一次膨脹過程就發生潰滅，VVF最大值為30%左右。與該處空泡潰滅的頻率相對應，在此處每一個振動周期內均出現脈沖壓力，其值約0.8 MPa。在試件邊緣處，空化效果相對較弱，VVF值較小，空泡在兩個振動周期內經過二次膨脹后潰滅，脈沖壓力約0.5 MPa。

從不同監測點的空泡體積變化曲線(圖6)和試樣振動速度變化曲線(圖3a)可看出，兩者具有相同的變化趨勢，表明試樣振動加速度決定了磁致伸縮儀空化流空化變化趨勢。當試樣的加速度為正值時，其方向為拉伸液體方向，試樣表面各處空泡體積組分不斷增加；反之，試樣表面空泡體積組分減小。

(a)中心

(b)徑向距離5 mm處

(c)邊緣圖5 試件表面不同位置壓力隨時間變化

(a)中心

(b)徑向距離5 mm處

(c)邊緣圖6 試件表面不同位置氣體體積組分隨時間變化

當t=0.775 ms時，靠近試樣的縱截面局部流場及試件表面的壓力和VVF分布云圖分別見圖7和圖8。當t=0.775 ms時，試件向下運動到平衡位置，此時加速度向下減小為0，速度達到最大。試件表面壓力和VVF圍繞試件中心近似呈環形分布，兩者分布相對應。此刻，隨壓力增大，靠近試件表面空泡發生集中潰滅，產生脈沖壓力，脈沖高壓分布在試件中心和間隔的環形區域。但在脈沖高壓區以外靠近試樣局部流場處仍存在較低負壓區，該區域有一定數量的空泡，VVF最大值可達33.5%。圖7和圖8表明即使在同一環形區域內，壓力和VVF也不是均勻分布的，存在無規律斷續脈動，這可能是由水流的脈動和水流運動影響空泡分布的均勻性及不均勻空泡潰滅所引起的。大于1 MPa的脈沖高壓區分布在圍繞試件中心半徑約1.2 mm的近似半球范圍內。由于試件表面壓力分布近似呈環狀分布，試樣空蝕形貌應呈現環狀分布的特征，這與大量的磁致伸縮儀空蝕實驗結果相一致。試件中心產生較高的脈沖壓力將使試件中心空蝕嚴重，但在外部的脈沖高壓區域，斷續出現的局部較高脈沖壓力將會引起不均勻分布的空蝕坑的出現。

圖9所示為相鄰兩周期內試樣向下運動到平衡位置時試樣表面沿徑向壓力分布。此時，試樣壓縮液體速度最大，試件表面空泡發生集中潰滅，產生脈沖壓力。從圖9可看出脈沖壓力按間隔環形區域分布，且隨試樣振動在相鄰的環形區域上交替出現。對于脈沖壓力交替出現的相鄰環形區域，其邊界有重疊。在環形區域重疊處，每一個振動周期均出現一次脈沖壓力，這與試樣中心和各環形區域內部每兩個振動周期出現一次脈沖壓力的頻次不同。從圖9可看出徑向距離為5 mm處恰好位于相鄰環形區域重疊處，圖5b中該處的壓力變化體現了這一特征。

圖7 t=0.775 s時縱向截面局部流場和試件表面的壓力分布云圖

圖8 t=0.775 ms時縱向截面局部流場和試件表面的VVF分布云圖

圖9 相鄰周期內空泡集中潰滅時試樣表面徑向壓力分布

圖10所示為相鄰兩個周期內試樣振動加速度向上過程中徑向VVF分布隨時間的變化。圖10中曲線1和曲線4分別為相鄰兩周期空泡集中潰滅時沿徑向VVF分布。曲線1→3和曲線4→6表明，對于相鄰兩周期試樣加速度指向拉伸液體的過程，靠近試樣局部流場的VVF值不斷增大。每一環形區域內空泡在兩個振動周期內經過二次膨脹后發生潰滅。通過對比曲線1和曲線4可看出試樣表面空泡潰滅發生在間隔環形區，且隨試樣振動，在相鄰環形區域上空泡發生交替潰滅，這與試樣表面脈沖壓力的分布特征相對應。在空泡潰滅區以外的間隔環形區域內存在一定數量的空泡，其體積被壓縮。

(a)

(b)圖10 相鄰周期內試件加速度向上過程中徑向VVF分布變化

圖11 t=0.775 ms時縱截面局部流場的壓力等值線

當t=0.775 ms時，縱截面局部流場的壓力等值線見圖11。圖11表明超聲壓力波不是以試件為中心呈球面波形傳播的。由于磁致伸縮儀空化流場的空化程度不均勻和空泡潰滅使試件表面附近局部流場形成高低壓區，試件表面可近似多個聲波發生源，各聲波傳播時相互疊加和干擾。

圖12所示為一個振動周期內在試件表面中心垂線上壓力分布的變化。變幅桿振動頻率為20 kHz，超聲波波長為70 mm。圖12表明在磁致伸縮儀超聲空化流場中，隨著距振動試樣距離增加，聲壓快速衰減，壓力的明顯波動只發生在距試件表面約20 mm內，約為波長的1/3，驗證了由變幅桿高頻振動產生的超聲波在空化流場傳播過程中能量快速衰減。

圖12 某一周期內變幅桿端面垂線上壓力變化

對于超聲空化流場中超聲波能量的快速衰減，除了超聲波的擴散衰減和由液體黏滯性引起的吸收衰減因素外，在靠近變幅桿端面附近的空化區域內，空泡的生成、膨脹等空化過程吸收超聲波能量及空泡對超聲波傳播的干擾等因素也會引起超聲波的能量衰減。

4 結論

(1) 通過數值模擬，驗證了磁致伸縮儀空化實驗試樣空蝕機理。變幅桿高頻振動引起試件表面附近局部流場發生空化。隨壓力波動，當試件表面附近空泡潰滅時產生脈沖壓力。在試樣中心區域空化程度和脈沖壓力作用效果顯著，最大脈沖壓力可達14 MPa。

(2)壓力和VVF在試樣表面近似呈環形分布。在同一環形區域內，壓力和VVF存在無規律斷續脈動。脈沖壓力在間隔的環形區域形成，且隨試樣振動在相鄰的環形區域上交替出現。

(3)在超聲空化流場中，試件表面可近似多個聲波發生源，各聲波傳播時相互疊加和干擾。聲壓隨距試樣距離增加快速衰減，壓力明顯波動只發生在距試件表面約20 mm內，約為波長的1/3。

[1] Philipp A, Lauterborn W. Cavitation Erosion by Single Laser-produced Bubbles[J]. Journal of Fluid Mechanics, 1998, 361: 75-116.

[2] 夏冬生, 張會臣, 于彥. 柴油機氣缸套冷卻水空化流的三維數值模擬[J]. 機械工程學報, 2011, 47(22):167-173. Xia Dongsheng, Zhang Huichen, Yu Yan. 3D Numerical Simulation of Cooling-water Cavitation Flow of Cylinder Liner for a Diesel Engine[J]. Journal of Mechanical Engineering, 2011, 47(22):167-173.

[3] Hattori S, Mikami N. Cavitation Erosion Resistance of Stellite Alloy Weld Overlays[J]. Wear, 2009, 267(11):1954-1960.

[4] Ahmed S M. Investigation of the Temperature Effects on Induced Impact Pressure[J]. Wear, 1998, 218(1):119-127.

[5] Benjamin T B, Ellis A T. The Collapse of Cavitation Bubbles and the Pressures Thereby Produced Against Solid Boundaries[J]. Philos. Trans. R. Soc., 1966, 260(1110):221-240.

[6] Leighton T G. The Acoustic Bubble[M]. London: Academic Press, 1997.

[7] M?rch K A. Erosion[M]. London: Academic Press, 1979.

[8] Bai Lixin, Xu Weilin, Zhang Faxing, et al. Cavitation Characteristics of Pit Structure in Ultrasonic Field[J]. Science in China (Series E: Technological Sciences), 2009,52(7):1974-1980.

[9] Burdin F, Tsochatzidis N A, Guiraud P, et al. Characterisation of the Acoustic Cavitation Cloud by Two Laser Techniques[J]. Ultrasonics Sonochemistry, 1999,6(1/2):43-51.

[10] Singhal A K, Athavale M M, Li H Y, et al. Mathematical Basis and Validation of Full Cavitation Model[J]. Journal of Fluid Engineering, 2002, 124(3):617-624.

[11] Menter F R. Two-equation Eddy-viscosity Turbulence Models for Engineering Applications[J]. AIAA Journal, 1994, 32(8):1598-1605.

(編輯 陳 勇)

3D Unsteady Numerical Simulation of Magnetostriction-induced Ultrasonic Cavitation Flow

Xia Dongsheng Sun Changguo Yu Yan Zhang Huichen

Dalian Maritime University, Dalia，Liaoning，116026

The magnetostriction-induced ultrasonic cavitation flow was numerically simulated by using a “Singhal full cavitation model” and SSTk-ωturbulence model with the dynamic mesh technique. The characteristics of ultrasonic cavitation flow and the cavitation erosion mechanism were revealed. The computational results show that pressure and vapor volume fraction (VVF) vary periodically on the specimen with vibration of the specimen. The lowest value of pressure fluctuation reaches the vapor pressure at the local flow field proximity to the specimen, where cavitation can occur. Pressure fluctuation presents a distinct pulse feature on the specimen surface, due to cavitation. Moreover, pressure and VVF were annularly distributed around the center on the specimen. However, they have irregular and discontinuous pulse in the annular zones. At the center of specimen, bubbles collapse leads to intense pulse pressure after twice oscillation, which can reach about 14 MPa. Pulse pressure was distributed in the interval annular zones and alternately occurs in the adjacent annular zones with vibration of the specimen. In the ultrasonic cavitation flow field, the specimen surface seems to be as several acoustic emission sources. The acoustic waves generate mutual superposition and interference. Pressure dramatically attenuates with increasing the distance to the specimen. The considerable pressure fluctuation occurs within the distance of about 20 mm to the specimen.

magnetostriction apparatus; computational fluid dynamics(CFD); cavitation; pulse pressure

2015-10-13

國家自然科學基金資助項目(51275064,50975036);遼寧省自然科學基金資助項目(2015020136)；中央高校基本科研業務費專項資金資助項目(3132013063)

TV131.3； O359.1

10.3969/j.issn.1004-132X.2016.22.014

夏冬生，男，1973年生。大連海事大學交通運輸裝備與海洋工程學院副教授。主要研究方向為機械摩擦磨損機理和潤滑、計算流體力學。發表論文20余篇。孫昌國，男，1979年生。大連海事大學交通運輸裝備與海洋工程學院講師。于 彥，女，1972年生。大連海事大學交通運輸裝備與海洋工程學院副教授。張會臣(通信作者)，男，1965年生。大連海事大學交通運輸裝備與海洋工程學院教授。