黃 瑩，高璞珍

（哈爾濱工程大學(xué) 核安全與仿真技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001）

在核能工程領(lǐng)域，兩相流動現(xiàn)象在核動力裝置中是經(jīng)常發(fā)生的，氣相在水中的存在形式以及運(yùn)動規(guī)律是影響兩相流動特性的關(guān)鍵所在，因此，氣泡運(yùn)動特性研究對涉及兩相運(yùn)動的核能工程領(lǐng)域具有重要的理論價(jià)值和實(shí)際應(yīng)用意義。

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

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

1 物理模型

1.1 幾何模型

計(jì)算域?yàn)橐粋€(gè)2D 的矩形計(jì)算區(qū)域，尺寸為40mm×80mm，在計(jì)算域中存在水平分布的兩個(gè)直徑為5mm、初始形狀為圓形的氣泡，網(wǎng)格為邊長0.1mm 的正方形網(wǎng)格。入口與出口邊界條件為速度入口和自由出流，本文模擬靜水情況，因此將入口速度設(shè)置為0。

1.2 搖擺模型

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

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

2 數(shù)學(xué)模型

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

2.1 VOF模型

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

氣體體積分?jǐn)?shù)方程為：

液相體積分?jǐn)?shù)方程與氣相體積分?jǐn)?shù)方程具有相同的形式，且兩相體積分?jǐn)?shù)滿足如下關(guān)系：

其中：αg和αl分別為氣相和液相的體積分?jǐn)?shù)；下標(biāo)l和g分別代表液相和氣相。不可壓縮流體的連續(xù)方程為：

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

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

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

F 為由附加的表面張力導(dǎo)致的動量源項(xiàng)，具有如下形式：

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

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

2.2 動網(wǎng)格模型

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

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

3 結(jié)果與討論

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

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

當(dāng)兩個(gè)氣泡一起運(yùn)動時(shí)，由于彼此的影響，其運(yùn)動特性較單氣泡情況有一定的區(qū)別，如當(dāng)兩氣泡距離較近時(shí)，其周圍的流場結(jié)構(gòu)會發(fā)生變化，其形狀及上升軌跡等都會隨之改變，且氣泡之間還可能發(fā)生碰撞、聚合等情況。

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

圖2為兩氣泡的變形過程及周圍的流場結(jié)構(gòu)?？捎^察到，兩氣泡在運(yùn)動過程中會使原本靜止的水產(chǎn)生流動，形成對稱分布的流場，通常表現(xiàn)為氣泡底部存在一個(gè)射流，在氣泡兩側(cè)均會有一個(gè)近似于圓形的流場結(jié)構(gòu)，且兩側(cè)流場的旋轉(zhuǎn)方向相反，氣泡左側(cè)周圍流體逆時(shí)針旋轉(zhuǎn)，右側(cè)周圍流體順時(shí)針旋轉(zhuǎn)。

在兩氣泡之間區(qū)域，兩氣泡產(chǎn)生的流場的水平分速度方向是相反的，速度疊加后，氣泡底部射流對氣泡的作用點(diǎn)會向疊加區(qū)域移動；氣泡位于此區(qū)域的部分受到的射流推力較大，上升速度較快，氣泡因兩側(cè)受力不平衡而發(fā)生旋轉(zhuǎn)，兩氣泡呈“八”字形分布（圖2d）。氣泡逐漸遠(yuǎn)離，彼此之間的影響變小，氣泡之間區(qū)域內(nèi)兩氣泡所產(chǎn)生的流場的水平分速度相互抵消效應(yīng)變小，射流對氣泡的作用點(diǎn)向氣泡中心移動（圖2e）。此時(shí)，兩氣泡靠近壁面?zhèn)葴?，處于水中壓力較大的位置，此處的空氣會向氣泡內(nèi)的低壓區(qū)移動，一段時(shí)間后，兩氣泡形狀及位置均處于近似對稱狀態(tài)（圖2f）。在氣泡位置由傾斜達(dá)到平衡的過程中，氣泡靠近壁面?zhèn)炔糠值乃俣容^快，使得氣泡兩側(cè)周圍水的速度有所提高，氣泡外側(cè)的速度大于位于疊加區(qū)域的部分，兩氣泡在相互靠近的同時(shí)會呈現(xiàn)反“八”字形分布（圖2g），氣泡的形狀及位置再次變得不平衡。此時(shí)氣泡間距離較近，中間區(qū)域流場疊加明顯，氣泡尾部射流對氣泡的作用點(diǎn)再次向疊加區(qū)域移動，氣泡位于疊加區(qū)域的部分因水流的推力作用，速度逐漸變大，氣泡形狀及位置再次達(dá)到平衡，上述過程構(gòu)成了一個(gè)位置變化周期。

表1 初始物性參數(shù)Table 1 Initial physical parameters

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

圖3為兩氣泡上升過程的速度分布情況，可看出，初始時(shí)，氣泡相互影響較弱，隨著形狀逐漸變得扁平，兩者之間距離變小，相互影響變強(qiáng)，如圖3b所示，兩氣泡之間區(qū)域流體下降速度大于另外兩側(cè)，從而產(chǎn)生一個(gè)低壓區(qū)，兩氣泡會向低壓區(qū)域靠近，兩氣泡相互吸引。當(dāng)兩氣泡呈“八”字形分布時(shí)，二者距離逐漸變大，如圖3d所示，兩氣泡流體速度變小，低壓區(qū)域消失，從而兩氣泡吸引作用變?nèi)酰瑲馀菰谏淞魍屏ο轮饾u分離。當(dāng)兩氣泡由“八”字形分布恢復(fù)到平衡位置時(shí)，如圖3f、g所示，氣泡間流體速度變大，氣泡相互吸引，兩氣泡再次靠近。

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

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

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

由于氣泡在上升過程中左右兩半部分受力不一致，使得周圍流體在氣泡擾動后，會出現(xiàn)漩渦，圖5為t=0.36s時(shí)管道內(nèi)的流場結(jié)構(gòu)。由圖4a可知，在0.36s時(shí)間內(nèi)，氣泡的位置關(guān)系經(jīng)歷了兩個(gè)周期，每個(gè)周期流體會產(chǎn)生兩對漩渦。

圖4 兩氣泡同時(shí)上升運(yùn)動位置周期Fig.4 Cyclicity of position relations between two rising bubbles

圖5 氣泡上升后產(chǎn)生的漩渦Fig.5 Vortex generated by rising bubble

3.2 搖擺對兩氣泡運(yùn)動規(guī)律的影響

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

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

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

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

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

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

4 結(jié)論

1）兩個(gè)水平并列分布的氣泡在上升過程會出現(xiàn)靠近-遠(yuǎn)離-再靠近-再遠(yuǎn)離的運(yùn)動現(xiàn)象，且二者的位置關(guān)系會呈現(xiàn)“八”字形與反“八”字形循環(huán)變化，模擬得到的結(jié)果與實(shí)驗(yàn)觀察到的現(xiàn)象符合良好。

2）兩氣泡并列上升時(shí)，每個(gè)氣泡都呈上下擺動狀上升，氣泡擺動使周圍水產(chǎn)生漩渦，在氣泡位置關(guān)系變化的一個(gè)周期內(nèi)，兩氣泡尾部會出現(xiàn)兩對漩渦。

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

［1］ LI Yong，ZHANG Jianping，F(xiàn)AN 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，F(xiàn)AN 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，F(xiàn)AN 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，F(xiàn)AN 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］ 高璞珍，龐鳳閣，王兆祥.核動力裝置一回路冷卻劑受海洋條件影響的數(shù)學(xué)模型［J］.哈爾濱工程大學(xué)學(xué)報(bào)，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高級應(yīng)用與實(shí)例分析［M］.北京：清華大學(xué)出版社，2010：32-35，140-182.

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