晃蕩條件下氣泡上升速度特性研究

宋禹林，譚思超，付學寬，李小輝

（1.哈爾濱工程大學核安全與仿真技術國防重點學科實驗室，黑龍江哈爾濱150001；2.中國艦船研究設計中心，湖北武漢430064）

晃蕩條件下氣泡上升速度特性研究

宋禹林1，譚思超1，付學寬1，李小輝2

（1.哈爾濱工程大學核安全與仿真技術國防重點學科實驗室，黑龍江哈爾濱150001；2.中國艦船研究設計中心，湖北武漢430064）

海浪、地震等條件都會導致核動力設備中儲存自由液面的設備產生劇烈的晃蕩。為了分析晃蕩條件下液體中的氣泡行為特性，本文運用CLSVOF（coupled level?set and volume?of?fluid）模型追蹤兩相流體交界面，數值模擬氣泡在液體中受余弦簡諧激勵和自由液面影響的上升過程。計算結果顯示余弦簡諧激勵降低了氣泡縱向運動速度，而橫向速度略大于激勵速度；自由液面的波動在縱向上會使氣泡上升速度驟然降低，甚至出現負速度，越接近自由液面橫向速度峰值越大。計算結果表明晃蕩運動的影響是余弦簡諧激勵與自由液面波動影響的疊加，晃蕩條件下自由液面波動是氣泡速度變化的主要因素，氣泡的本身屬性起次要作用。

上升氣泡；晃蕩；CLSVOF；大空間

行駛在海洋中的船舶、運行在陸地上的核電站，在遇到風浪、地震等外界激勵時，會使存在自由液面的容器產生劇烈的晃蕩。強烈的晃蕩產生的復雜流場必然對氣泡的運動有著不可忽視的影響。正常工況下運行的蒸汽發生器、除氧器還有事故工況下未停止運行的反應堆都存在大量的氣泡，氣泡的大小、形狀及運動規律等對氣液傳熱傳質都有這種影響，因此晃蕩對氣泡運動行為的影響或危機設備的安全運行。而這方面的研究卻未見諸公開報道，因此有必要對晃蕩條件下氣泡的運動行為特性進行研究。

1 數值計算方法與靜水中氣泡的上升

氣泡運動的數值模擬主要通過界面追蹤法來實現，Level Set和VOF都是優秀的界面跟蹤方法，C.W.Hirt［1］首次將VOF方法應用到數值模擬晃蕩現象中，Sussman M［2］等將Level Set方法應用于數值模擬氣泡的上升運動。Son G［3］等完善了CLSVOF方法，相比其他兩相流界面追蹤方法，CLSVOF方法計算得更精準并且收斂的更快。耦合Level Set方法與VOF方法關鍵的過程在于相函數的耦合方法。相函數耦合又分為以Level Set方法為主還是以VOF方法為主，這里選以Level Set為主的Level Set-相界面法為例，式（1）是一種最基本的耦合方法，利用Heaviside函數建立了Level Set函數φ（x，y，z，t）和VOF流體體積函數f（Ω，t）的關系［4］：

運動速度主要取決于氣泡的尺寸和形狀，在氣泡的上升過程中，由于內外壓差以及速度差氣泡體積會被壓縮，根據Eotvos數的不同會形成不同的形狀，依照文獻［5］的實驗結論，直徑為0.02 m的氣泡在上升達到最終速度時應該為扁橢球，此時可以用Wairegi與Grace在實驗中擬合的經驗公式［6］來計算氣泡上升速度。上升氣泡在靜水中速度隨時間變化的計算結果如圖1所示，最終速度為0.1884 m／s，與文獻［6］的實驗結果相符合，氣泡x方向的速度滿足氣泡螺旋式上升這一普遍現象。

圖1 靜水中氣泡上升速度Fig.1 Rising velocity of bubbles in static water

2 晃蕩對大空間內液體流場的影響

晃蕩條件下氣泡運動速度的變化，最直接的影響因素是氣泡周圍的流場。而晃蕩本質上是由外界激勵驅動的，因此晃蕩條件下容器內的流場主要受到外界激勵和自由液面波動這2方面的綜合影響。本章主要研究外界激勵和自由液面對容器內流場的影響。

2.1 余弦激勵條件下的流場

首先研究寬0.4 m、高0.6 m充滿液體的容器在受到余弦簡諧激勵時的流場變化。觀察1／2周期內流場x方向速度，如圖2所示，流體受到簡諧激勵的作用，容器上下壁面附近流體跟隨激勵運動，流場內速度由容器上下壁面向流場中心傳遞。

容器在受簡諧激勵運動的過程中，流場狀態在大部分時間都處于圖2所示的狀態。從云圖中可以看出流場x方向速度變化不大，另外如圖3所示x方向速度的標準差，最大值僅為0.009 6。

圖2 流場x方向速度云圖Fig.2 Velocity contour diagram of direction x in the flow field

圖3 流場x方向速度標準差Fig.3 Velocity standard deviation of direction x in the flow field

觀察流場y方向速度，流體存在y方向速度是因為流體本身具有可壓縮性，而水的可壓縮性很低，因此y方向速度不會太大，如圖4所示為y方向速度占x方向速度的百分比，最大值僅為0.000 5%，因此表征流場速度時y方向速度可以忽略，僅用x方向平均速度表示即可。而流場x方向速度又可以近似為簡諧激勵速度，因此氣泡x方向運動速度直接與簡諧激勵速度相關。

圖4 流場y方向速度占x方向速度百分比Fig.4 Ratio of velocity in direction y to direction x in the flow field

2.2 自由液面波動條件下的流場

大空間內無自由液面的流體受到簡諧激勵時，流場x方向速度與激勵速度相同，并且流場內速度規則分布。存在自由液面的流場在受到簡諧激勵時，流場內速度分布與無自由液面本質相同，均是由上下容器壁面向中心傳遞。此時中心位置不再是幾何中心，而是接近自由液面的位置。如圖5所示，自由液面將流場分為2部分，上部分氣空間流場流速明顯大于下部分液空間流場流速，下部液體空間流場流速略大于0.112 2 m／s的激勵速度。

圖5 流場x方向速度云圖Fig.5 Velocity contour diagram of x direction in the flow field

如圖6所示隨著速度逐漸向流場中心傳遞，中心速度也會大于上下容器壁面附近流場速度。由于自由液面的存在使得越接近自由液面，流場內x方向的運動速度就會越大并且會超過激勵速度。

圖6 流場x方向速度云圖Fig.6 Velocity contour diagram of direction x in the flow field

無自由液面的流體受到簡諧激勵時流場x方向速度均方差非常小，可以用平均速度表征整個流場內的速度。但存在自由液面時，流場內x方向速度均方差較大，如圖7所示均方差最低也大于0.020 5，是無自由液面時最大值0.009 6的2.14倍，最大值0.176 2是無自由液面最大值的18.35倍。均方差較大，表示數據偏移平均值程度較大，此時無法用流場x方向速度平均值表征流場內流體x方向速度。

圖7 流場x方向速度標準差Fig.7 Velocity standard deviation of direction x in the flow field

無自由液面的流體在受到簡諧激勵時，y方向的速度主要由液體的可壓縮性引起，水的可壓縮性較小導致y方向速度也很小，在表征流場流速時可以忽略y方向速度。但大空間內有自由液面的流體y方向速度很大，受到簡諧激勵5 s內最大速度能夠達到0.49 m／s，最小速度也達到了-0.52 m／s。如圖8所示為0.2 s是的流場y方向速度云圖。云圖中的0速度線位置為流場中線，流場被0速度線分為2部分，一部分速度為正，另一部分為負，2部分的中心位置皆為容器壁面約0.4 m處，略高于自由液面初始高度，無論自由液面如何波動，中心位置都只在極小范圍內變化。

圖8 流場x方向速度云圖Fig.8 Velocity contour diagram of direction x in the flow field

如圖9所示，整1／2周期時右側液面達到最高位置，此時激勵速度為0，但流場內y方向最大速度為0.29 m／s、最小速度-0.31 m／s，0速度分界線依然在流場中線處，2部分流場的中心位置與0.2 s相比偏移量極小。由此可知隨著氣泡的上升受到自由液面波動引起的y方向沖擊越大，因此越接近自由液面氣泡的y方向速度波動越大。流場y方向的0速度線近似為垂直中線，偏離0速度線越遠，受y方向速度影響也越明顯。隨著氣泡的上升，氣泡x方向速度波動越大，偏離0速度線越遠，距離水面越近，流場y方向上速度對氣泡的影響逐漸加劇，并且靜水中氣泡上升的最終速度為0.188 4 m／s，單y方向速度最小-0.31 m／s，因此氣泡上升速度會出現負值。

圖9 流場y方向速度云圖Fig.9 Velocity contour diagram of direction y in the flow field

3 晃蕩對上升氣泡速度的影響

3.1 余弦簡諧激勵對氣泡上升速度的影響

首先考慮單純受到外界激勵的情況（即無自由液面），對受到式（2）的余弦簡諧激勵下氣泡的上升運動進行數值模擬，容器寬0.4 m、高0.6 m、液位0.6 m，氣泡初始直徑0.02 m，距離容器底部中心位置0.02 m、初始速度為0，振幅取六分之底邊長A＝0.1 m，周期選取1.32 s。

數值模擬結果如圖10、圖11所示，由于整個系統中慣性力起主導作用，氣泡x方向速度與簡諧激勵速度相比，主要差別在于氣泡運動速度峰值略大于激勵速度，氣泡運動速度變化先于激勵速度變化。由于慣性力的加入氣泡y方向受到的阻力加大，上升氣泡y方向速度小于靜水中氣泡的上升速度，并且難以達到動平衡狀態。

液體流速變化與激勵速度同步，氣泡速度變化先于激勵速度。這是由于空氣和水的密度差也為氣泡提供驅動力，使氣泡運動的驅動力大于激勵驅動，密度差越大驅動力越大，密度差越小氣泡運動速度越接近流體，隨著氣體密度的增加差距逐漸減小，若與液體密度相同則氣泡運動也與激勵速度完全同步，如圖12所示。

圖10 氣泡y方向速度Fig.10 Velocity in direction y of bubbles

圖11 氣泡x方向速度Fig.11 Velocity in direction x of bubbles

圖12 不同密度的氣泡速度Fig.12 Velocities of bubbles with different densities

3.2 自由液面波動對氣泡上升速度的影響

自由液面在外界激勵的影響下會產生強烈的波動，為研究大空間內自由液面波動對氣泡上升速度影響，現在設定自由液面高度為0.4 m，其他條件不變，上升氣泡速度如圖13所示。

氣泡y方向運動速度不再是逐漸增加到最終速度保持穩定，而是先逐漸增加，然后因此受到自由液面波動的沖擊阻礙氣泡上升，導致速度下降，沖擊減弱后速度再次增加，而后不斷重復這一過程，越接近自由液面速度下降時越劇烈，甚至產生負速度。氣泡x方向運動速度的變化規律總體上與附加慣性力條件下相似，但受到自由液面波動的疊加，越接近自由液面氣泡擺動速度越大，一周期后氣泡擺動速度始終大于激勵速度。

圖13 晃蕩條件下氣泡的運動速度Fig.13 Velocity of bubbles under sloshing condition

3.3 氣泡速度影響因素的頻譜分析

本節將利用頻譜分析的手段來驗證晃蕩對氣泡運動的影響，首先對有自由液面的流場進行分析。由無激勵條件下的氣泡運動速度可以預估晃蕩條件下的氣泡速度，將氣泡抽象成質點，該質點在預估速度下上升，對該質點位置上流場速度的時間序列進行頻譜分析，該速度的振幅頻譜如圖14所示，本算例中流場速度的影響因素只能是外界激勵和自由液面波動，0.76 Hz峰值對應外界激勵，0.2 Hz峰值對應自由液面波動。因此可以說在本文中簡諧激勵條件下自由液面的波動，對直徑為0.02 m的氣泡上升過程產生的影響頻率為0.2 Hz，并且越接近自由液面0.2 Hz頻譜峰值越大。

圖14 自由液面波動的頻率Fig.14 Frequency of the free surface fluctuation

對靜水中氣泡上升速度進行傅里葉變換，靜水中氣泡速度的振幅頻譜如圖15所示。2個大峰值的位置在1～2.5 Hz，說明氣泡所受浮力、升力以及氣泡周圍擾流等由氣泡自身特性引起的外力對氣泡x方向運動速度的影響，主要集中在1～2.5 Hz這個范圍內。

圖15 靜水中氣泡速度的振幅頻譜Fig.15 Amplitude frequency spectrum of bubble veloci?ty in static water

圖16 為晃蕩條件下氣泡速度的振幅頻譜圖。圖中峰值B對應頻率0.76 Hz，表示余弦簡諧激勵，即對應附加慣性力對氣泡運動速度的影響；又根據先前的討論，峰值C、D對應氣泡所受阻力、浮力等因素的影響。峰值A對應自由液面波動對氣泡運動速度的影響。

圖16 晃蕩條件下氣泡速度的振幅頻譜Fig.16 Amplitude frequency spectrum of bubble veloci?ty under sloshing condition

綜上所述，自由液面波動劇烈時對氣泡運動速度的影響強于附加慣性力產生的影響，氣泡自身特性對氣泡運動速度的影響最弱，自由液面波動影響最強烈。

4 結論

本文采用CLSVOF方法追蹤兩相交界面，數值模擬晃蕩條件下的氣泡上升運動。分別得出余弦簡諧激勵和自由液面波動影響下上升氣泡的速度變化規律，然后分析了余弦簡諧激勵和自由液面波動影響下的流場變化，從流場變化的角度分析氣泡速度的變化。最后從頻譜分析的角度分析氣泡上升過程中受到的影響，分析出了晃蕩對氣泡運動速度的影響機理，總結出如下結論：

1）晃蕩對氣泡運動速度變化的影響分為自由液面波動產生的影響、附加慣性力產生的影響以及氣泡自身特性的影響，其中自由液面波動影響最大，氣泡自身特性影響最小。

2）余弦簡諧激勵影響下流場y方向速度極小可以忽略，x方向速度平均值可以表征流場流動速度，并且與激勵速度相等。氣泡x方向速度大于激勵速度，并且變化先于激勵速度，差距隨著氣泡密度的增加而減小。氣泡y方向速度低于靜水中速度，并且難以達到動平衡狀態。

3）自由液面的波動使流場獲得較大的y方向速度，x方向速度平均值的峰值也大于簡諧激勵。隨著氣泡的上升x方向速度峰值逐漸增大，y方向速度在收到沖擊時會下降，越接近自由液面速度下降時越劇烈，甚至產生負速度。

4）CLSVOF法雖然計算準確但計算量大消耗時間，在未來的研究中可以引入邊界元法或通過實驗進一步深入研究。

［1］DODGE F T，TEXAS S A.The new dynamic behavior of liq?uids in moving containers［R］／／San Antonio：Southwest Re?search Institute，2000.

［2］SUSSMAN M，SMEREKA P，OSHER S.A level set method for computing solutions to incompressible two?phase flows［J］.Journal of Computational Physics，1994，114（1）：146?159.

［3］SON G，HUR N.A coupled level set and volume?of?fluid method for the buoyancy?driven motion of fluid particles［J］.Numerical Heat Transfer，Part B：Fundamentals：An Inter?national Journal of Computation and Methodology，2002，42（6）：523?542.

［4］CHAKRABORTY I，BISWAS G，GHOSHDASTIDAR P S.A coupled level?set and volume?of?fluid method for the buoy?ant rise of gas bubbles in liquids［J］.International Journal of Heat and Mass Transfer，2013，58：240?259.

［5］BHAGA D，WEBER M E.Bubbles in viscous liquids：shapes，wakes and velocities［J］.Journal of Fluid Mechan?ics，1981，105：61?85.

［6］WAIREGI T，GRACE J R.The behavior of large drops in im?miscible liquids［J］.International Journal of Multiphase Flow，1976，3（1）：67?77.

［7］TADASHI W.On the numerical approach for simulating re?actor thermal hydraulics under seismic conditions［J］.An?nals of Nuclear Energy，2012，49：200?206.

［8］劉云龍，張阿漫，王詩平，等.基于邊界元法的氣泡同波浪相互作用研究［J］.物理學報，2012，61（22）：269?276.LIU Yunlong，ZHANG Aman，WANG Shiping，et al.Re?search on interaction between bubble and surface waves based on BEM［J］.Acta Physica Sinica，2012，61（22）：269?276.

［9］倪寶玉，張阿漫，張忠宇，等.復雜流場下氣泡界面不穩定現象的數值模擬方法綜述［J］.船舶力學，2012，4：470?483.NI Baoyu，ZHANG Aman，ZHANG Zhongyu，et al.Review on the numerical simulation methods of the instability of bub?ble interface under complex flow field［J］.Journal of Ship Mechanics，2012，18（4）：470?483.

Study on the characteristics of single bubble rising velocity under sloshing conditions

SONG Yulin1，TAN Sichao1，FU Xuekuan1，LI Xiaohui2

（1.Fundamental Science on Nuclear Safety and Simulation Technology Laboratory，Harbin Engineering University，Harbin Engineering University，Harbin 150001，China；2.China Ship dEvelopment and Design Center，Wuhan 430064，China）

Ocean waves and seismic waves can cause violently sloshing of the free liquid surface in a container.In order to analyze the behavior of rising bubble in liquid under the condition of sloshing，this paper uses the CLSVOF（coupled level?set the volume?and?fluid）model to simulate the two?phase fluid interface，obtaining the rising bub?ble in the liquid with cosine harmonic excitation condition or wave of free surface condition.Calculation results show that cosine harmonic excitation reduces the bubble longitudinal velocity transverse velocity is slightly bigger than the motivation speed and the wave of free surface will make the bubble rising velocity plummet and even lead to nega?tive speed.The lateral velocity peak increases gradually closer to the free surface.Calculation results prove that sloshing is the outcome of superposition with cosine harmonic excitation and wave of free surface and the density difference between air and water is one of the driving forces of bubble motion.In sloshing condition，fluctuation of free surface is the main factor of changes in bubble velocity，which proved that the attribute of bubble is the wea?kest factor.

rising bubble；sloshing；CLSVOF；large space

10.3969／j.issn.1006?7043.201404023

TL334

：A

：1006?7043（2015）06?0865?06

http：／／www.cnki.net／kcms／detail／23.1390.u.20150428.1115.017.html

2014?04?07.網絡出版時間：2015?04?28.

黑龍江省青年學術骨干支持計劃資助項目（1254G017）；核安全與仿真技術國防重點學科實驗室基金資助項目（HEUFN1305）.

譚思超（1979?），男，教授，博士生導師.

譚思超，E?mail：tansichao＠hrbeu.edu.cn.