水平分布兩氣泡之間相互作用及搖擺對其影響的數值模擬

黃 瑩，高璞珍

（哈爾濱工程大學 核安全與仿真技術國防重點學科實驗室，黑龍江 哈爾濱 150001）

在核能工程領域，兩相流動現象在核動力裝置中是經常發生的，氣相在水中的存在形式以及運動規律是影響兩相流動特性的關鍵所在，因此，氣泡運動特性研究對涉及兩相運動的核能工程領域具有重要的理論價值和實際應用意義。

近年來，隨著計算機技術和數值方法的發展，數值計算越來越廣泛地應用到氣泡動力學的研究，文獻［1-4］應用VOF 方法模擬了高壓下單個氣泡上升過程中的變形情況。Chen等［5］應用VOF方法模擬了單個氣泡在靜止水中的上升及其與氣水交界面的相互作用過程。

在現有的氣泡運動數值模擬中，多數是針對單個氣泡運動規律的研究，對于水平分布的兩個氣泡間相互作用及氣泡對的上升過程的模擬則較少，VOF方法在保證具有較高精度的前提下，具有較小的計算量，相較其他數值方法更易實現，其在氣泡研究領域具有不可替代的作用。本文擬采用VOF 中的PLIC 界面捕捉方法，在動量控制方程中加入表面張力項，對重力作用下兩個氣泡的上升過程及相互作用規律進行二維數值模擬研究。

1 物理模型

1.1 幾何模型

計算域為一個2D 的矩形計算區域，尺寸為40mm×80mm，在計算域中存在水平分布的兩個直徑為5mm、初始形狀為圓形的氣泡，網格為邊長0.1mm 的正方形網格。入口與出口邊界條件為速度入口和自由出流，本文模擬靜水情況，因此將入口速度設置為0。

1.2 搖擺模型

由于實際的海洋條件較復雜，一般可繞著管道的初始平衡位置發生6 個自由度的振蕩運動，包括前后搖擺、左右搖擺、上下起伏，同時由于搖擺，還會出現管道傾斜等現象。本文在進行搖擺模擬時，做了適當的簡化，將搖擺模型簡化為管道繞軸線做角度變化為θ=θmsin（2πt／T）（其中θm=15°，為最大擺角幅值；T=6s，為搖擺周期）的簡諧搖擺，如圖1所示，搖擺中心位于管道左下角，并做以下假設：平衡位置為穩態時的位置；恢復力由波浪提供；搖擺時只發生角位移。此假設下，管道的位移角、角速度、角加速度均呈正弦規律變化［6］。

圖1 搖擺示意圖Fig.1 Swing diagram

2 數學模型

本文基于Fluent平臺，采用VOF 模型及PLIC界面捕捉方法，結合考慮了表面張力的動量方程，對靜水及搖擺情況下的氣泡間相互作用特性進行數值模擬，其基本數學模型構想如下。

2.1 VOF模型

VOF方法的基本原理是通過研究網格單元中流體和網格體積比函數來確定自由界面，追蹤流體的變化，而非追蹤自由界面上質點的運動，在每個控制體積內，兩相的體積分數之和為1。在單元中，如果第n相流體的體積分數記為αn，那么就會出現下面3個可能的情況：αn=0，第n相流體在單元中是空的；αn=1，第n相流體在單元中是充滿的；0＜αn＜1，單元中包含了第n 相流體和其他1相或多相流體的界面［7］。

氣體體積分數方程為：

液相體積分數方程與氣相體積分數方程具有相同的形式，且兩相體積分數滿足如下關系：

其中：αg和αl分別為氣相和液相的體積分數；下標l和g分別代表液相和氣相。不可壓縮流體的連續方程為：

其中，u為流體速度，m／s。

考慮表面張力的動量方程為：

其中：ρ為兩相平均密度，kg／m3；p 為壓力，N；μ 為兩相平均動力黏度，Pa·s；D 為應力張量，滿足如下關系：

F 為由附加的表面張力導致的動量源項，具有如下形式：

其中：κ為界面曲率；σ為表面張力系數，N／m。

出現在輸運方程中的屬性是由存在于每一控制體體積中的分相決定的，即：

2.2 動網格模型

本文通過UDF 編程手段，在動網格模型下實現搖擺工況的模擬，在計算過程中，除了要考慮流體自身的運動，還需附加上邊界的運動，這使得其控制方程會與其他情況有一定的區別。動網格模型中，在任一控制體V 內，其邊界是運動的，動量方程［8］為：

其中：vx、vy分別為控制體沿x 軸及y 軸方向的運動速度，m／s；a 為 控制體加速度，m／s2；f 為作用在單位質量流體上的質量力，N。

3 結果與討論

3.1 上升過程氣泡之間的相互作用

在一個40 mm×100 mm 的通道內（介質為水），距容器底部5mm 處水平并排分布兩個直徑為5mm 的球形氣泡，兩氣泡中心的初始水平距離為5mm。模擬中液體與氣體的物性參數列于表1。

當兩個氣泡一起運動時，由于彼此的影響，其運動特性較單氣泡情況有一定的區別，如當兩氣泡距離較近時，其周圍的流場結構會發生變化，其形狀及上升軌跡等都會隨之改變，且氣泡之間還可能發生碰撞、聚合等情況。

氣泡對水的擾動有一定的影響范圍，只有氣泡位于另一氣泡的影響范圍之內，才會產生相互作用。在作用過程中發現，氣泡首先會有一相互吸引的趨勢，然后又會發生分離，過程中氣泡之間存在一最小距離，若初始兩氣泡的間距很小，即小于此最小距離時，氣泡會發生接觸，進行融合。本文只針對水平分布的兩個氣泡在靜水中上升過程展現的運動特性進行研究，不考慮碰撞及聚合過程。

圖2為兩氣泡的變形過程及周圍的流場結構。可觀察到，兩氣泡在運動過程中會使原本靜止的水產生流動，形成對稱分布的流場，通常表現為氣泡底部存在一個射流，在氣泡兩側均會有一個近似于圓形的流場結構，且兩側流場的旋轉方向相反，氣泡左側周圍流體逆時針旋轉，右側周圍流體順時針旋轉。

在兩氣泡之間區域，兩氣泡產生的流場的水平分速度方向是相反的，速度疊加后，氣泡底部射流對氣泡的作用點會向疊加區域移動；氣泡位于此區域的部分受到的射流推力較大，上升速度較快，氣泡因兩側受力不平衡而發生旋轉，兩氣泡呈“八”字形分布（圖2d）。氣泡逐漸遠離，彼此之間的影響變小，氣泡之間區域內兩氣泡所產生的流場的水平分速度相互抵消效應變小，射流對氣泡的作用點向氣泡中心移動（圖2e）。此時，兩氣泡靠近壁面側滯后，處于水中壓力較大的位置，此處的空氣會向氣泡內的低壓區移動，一段時間后，兩氣泡形狀及位置均處于近似對稱狀態（圖2f）。在氣泡位置由傾斜達到平衡的過程中，氣泡靠近壁面側部分的速度較快，使得氣泡兩側周圍水的速度有所提高，氣泡外側的速度大于位于疊加區域的部分，兩氣泡在相互靠近的同時會呈現反“八”字形分布（圖2g），氣泡的形狀及位置再次變得不平衡。此時氣泡間距離較近，中間區域流場疊加明顯，氣泡尾部射流對氣泡的作用點再次向疊加區域移動，氣泡位于疊加區域的部分因水流的推力作用，速度逐漸變大，氣泡形狀及位置再次達到平衡，上述過程構成了一個位置變化周期。

表1 初始物性參數Table 1 Initial physical parameters

圖2 兩氣泡周圍流線分布及變形過程Fig.2 Streamline shape and deformation process of two rising bubbles

圖3為兩氣泡上升過程的速度分布情況，可看出，初始時，氣泡相互影響較弱，隨著形狀逐漸變得扁平，兩者之間距離變小，相互影響變強，如圖3b所示，兩氣泡之間區域流體下降速度大于另外兩側，從而產生一個低壓區，兩氣泡會向低壓區域靠近，兩氣泡相互吸引。當兩氣泡呈“八”字形分布時，二者距離逐漸變大，如圖3d所示，兩氣泡流體速度變小，低壓區域消失，從而兩氣泡吸引作用變弱，氣泡在射流推力下逐漸分離。當兩氣泡由“八”字形分布恢復到平衡位置時，如圖3f、g所示，氣泡間流體速度變大，氣泡相互吸引，兩氣泡再次靠近。

圖3 兩氣泡上升過程速度分布Fig.3 Velocity distribution of two rising bubbles

上升過程中，兩氣泡會出現靠近-分離-再靠近-再分離的運動現象，且過程中伴隨著兩氣泡呈現“八”字形與反“八”字形循環的現象。圖4a為模擬所得兩氣泡位置關系的變化過程，其中8組數據從下至上依次為t=0.04、0.08、0.12、0.16、0.22、0.28、0.32、0.36s時刻氣泡的位置。

圖4b為實驗觀察到的氣泡運動位置變化，圖4c為實驗段示意圖。雖然實驗中氣泡沿水平方向運動，但在相同的水平高度下，兩氣泡之間的作用可不考慮重力，可體現出兩并列氣泡沿同一方向運動時，其位置關系的變化規律。從圖中可看出，模擬結果與實驗觀察到的結果一致。

由于氣泡在上升過程中左右兩半部分受力不一致，使得周圍流體在氣泡擾動后，會出現漩渦，圖5為t=0.36s時管道內的流場結構。由圖4a可知，在0.36s時間內，氣泡的位置關系經歷了兩個周期，每個周期流體會產生兩對漩渦。

圖4 兩氣泡同時上升運動位置周期Fig.4 Cyclicity of position relations between two rising bubbles

圖5 氣泡上升后產生的漩渦Fig.5 Vortex generated by rising bubble

3.2 搖擺對兩氣泡運動規律的影響

在搖擺工況下，氣泡運動的獨特性是由兩方面因素決定的：1）氣泡受到搖擺運動產生的附加慣性力的作用；2）流體隨壁面搖擺運動產生的速度對氣泡的影響，氣泡在水中的變形過程主要受水的阻力影響。氣泡間的相互作用也是通過水為介質實現的，因此，在一般情況下，后一因素的作用會明顯大于前一因素的作用。

由于搖擺中心位于管道左下角，初始時，搖擺會使水具有向左上方流動的速度，從而影響氣泡上升所引起的對流，使氣泡周圍的流場與非搖擺情況不同。搖擺會改變管道內流體原有的速度場及壓力場結構，使二者均具有不對稱的結構，氣泡底部原本向上流動的射流在疊加了搖擺所產生的速度后，其射流方向會向通道左上方偏移，如圖6所示。

圖6 搖擺情況下氣泡周圍流線形狀Fig.6 Streamline shape around two rising bubbles in swinging condition

圖7為搖擺條件下氣泡的變形過程，與非搖擺條件氣泡的變形過程（圖2）對比可看出，搖擺工況下，兩氣泡在上升過程中會變得更加扁平，這是因為氣泡兩側的速度旋轉區域內的流體受搖擺向上速度分量的影響，旋轉區域會向上移動，從氣泡兩側下方移至氣泡兩側，流體的旋轉會產生一個低壓區域，這就意味著較非搖擺情況，氣泡更易向兩側低壓區域伸展，形狀更加扁平。

從圖7可發現，搖擺情況下兩個水平分布的氣泡由于形狀變得更加扁平，兩者距離會逐漸變小，所能達到的最小距離要小于非搖擺條件下的值，更容易相互影響，增大了聚合的幾率，擴大了可發生聚合條件的初始距離值的范圍。但在達到最小距離后，由于兩者距離搖擺中心的距離不同，受搖擺影響的效果也就不同，兩氣泡的速度大小及方向都會有所區別，使得兩氣泡逐漸分開，相互作用減弱，兩氣泡在上升過程所展現出的周期性變化會逐漸消失。但在搖擺運動初期，兩氣泡仍可展現出周期性的位置關系變化過程。

圖7 搖擺情況下水平分布兩氣泡上升的變形過程Fig.7 Deformation process of two rising bubbles in swinging condition

4 結論

1）兩個水平并列分布的氣泡在上升過程會出現靠近-遠離-再靠近-再遠離的運動現象，且二者的位置關系會呈現“八”字形與反“八”字形循環變化，模擬得到的結果與實驗觀察到的現象符合良好。

2）兩氣泡并列上升時，每個氣泡都呈上下擺動狀上升，氣泡擺動使周圍水產生漩渦，在氣泡位置關系變化的一個周期內，兩氣泡尾部會出現兩對漩渦。

3）搖擺會使氣泡周圍的流場結構變得不對稱，氣泡上升過程中變形過程更加明顯，氣泡變得更加扁平，兩氣泡所能達到的最小距離變小，但兩氣泡受搖擺影響不同，兩氣泡會逐漸分開，搖擺會破壞兩氣泡上升過程位置關系所展現出的周期性。

［1］ LI Yong，ZHANG Jianping，FAN Liangshi.Discrete-phase simulation of single bubble rise behavior at elevated pressures in a bubble column［J］.Chemical Engineering Science，2000，55（20）：4 597-4 609.

［2］ LI Y，ZHANG J P，FAN L S.Numerical studies of bubble dynamics in gas-liquid-solid fluidization at high pressures［J］.Power Technology，2001，116（2）：246-260.

［3］ YANG G Q，DU Bing，FAN L S.Bubble formation and dynamics in gas-liquid-solid fluidization：A review［J］.Chemical Engineering Science，2007，62（1）：2-27.

［4］ ZHANG Jianping，LI Yong，FAN Liangshi.Numerical studies of bubble and particle dynamics in a three-phase fluidized bed at elevated pressures［J］.Powder Technology，2000，112（1）：46-56.

［5］ CHEN L，LI Y.A numerical method for twophase flows with an interface［J］.Environmental Model and Software，1998，13（3-4）：247-255.

［6］ 高璞珍，龐鳳閣，王兆祥.核動力裝置一回路冷卻劑受海洋條件影響的數學模型［J］.哈爾濱工程大學學報，1997，18（1）：24-27.GAO Puzhen，PANG Fengge，WANG Zhaoxiang.Mathematical model of primary coolant in nuclear power plant influenced by ocean conditions［J］.Journal of Harbin Engineering University，1997，18（1）：24-27（in Chinese）.

［7］ 江帆，黃鵬.Fluent高級應用與實例分析［M］.北京：清華大學出版社，2010：32-35，140-182.

［8］ 孔瓏.工程流體力學［M］.3版.北京：中國電力出版社，2008：225-230.