船海結構物水氣泡混合流機理及數值模擬方法研究進展

張曉嵩，和康健，萬德成*,3

1 上海交通大學 船海計算水動力學研究中心，上海 200240

2 上海交通大學 船舶海洋與建筑工程學院，上海 200240

3 浙江大學 海洋學院，浙江 舟山 316021

0 引 言

在海洋環境中，船舶與海洋工程結構物在與周圍流體相互作用時，會出現明顯的白色泡沫流動，即為船海結構物水氣泡混合流。在實際作業過程中，船海結構物水氣泡混合流廣泛存在且形式多樣。對于運動中的船舶，由于船體與流體之間的強烈撞擊，船艏部會出現興波破碎，產生氣體卷吸和液體飛濺，從而形成水氣泡混合流[1]；在船艉部，由于船艉后方的流動分離，以及螺旋槳空泡流動與湍流的相互作用，會引起劇烈的水氣泡多相流的相互滲透與混合，從而形成大范圍擴散分布的白色氣泡尾跡[2]。船舶周圍產生的水氣泡混合流主要由船舶的航行運動引起，而海洋平臺周圍出現的水氣泡混合流則是由平臺與波浪間劇烈的相互作用引起。例如，波浪砰擊平臺立柱或平臺下甲板時，會導致水氣相互摻混與液滴飛濺霧化，形成水氣泡混合流[3]；在極端海況作用下，波浪沖擊平臺主體會產生甲板上浪和劇烈的沖刷流動，從而引起劇烈的水氣摻混，形成水氣泡混合流[4]。以上現象均為船海結構物水氣泡混合流的表現形式。作為船舶與海洋結構物等眾多自由液面工程實際問題的重要組成部分，這種水氣泡混合流是普遍存在、易于觀察的，但也難以系統描述和分析。

從流動組成上看，水氣泡混合流現象表現為：在水面以下存在大量的氣體卷吸進而演化成多尺度的氣泡，水面以上則存在大范圍的液體飛濺進而演化成霧化液滴，而氣泡和液滴同時在大變形自由面附近發生兩相劇烈的摻混，形成水、氣、泡、霧耦合的流動狀態。

水氣泡混合流的生成演化包含著劇烈的多相流動和多物理場的力學特性，對船舶結構物的水動力學性能、極端載荷、流激噪聲、尾跡形態等有著重要影響，重點體現在以下方面：

1） 水氣泡混合流是船海結構物直接與周圍流體劇烈地相互作用而生成，會引起船海結構物能量損失或是載荷變化，直接導致船舶快速性能和操縱性能的降低，或海洋平臺結構物耐波性能和穩定性能的改變[5]。

2） 水氣泡混合流中的分散多相流動與湍流相互作用，會引發流動劇烈脈動，改變船海結構物周圍流激噪聲的輻射特性，并會對聲學探測設備的信號傳輸產生影響[6]。

3） 水氣泡混合流會導致波浪砰擊結構物的載荷特性與搖蕩運動特性發生變化，飛濺高度增加，載荷擴散范圍增大，并加劇甲板上浪現象[7]。

4） 艦船尾跡中的水氣泡混合流容易被探測和追蹤，致使艦船隱身性能大大降低[8]。

由此可見，船海結構物水氣泡混合流與船海工程裝備的綠色節能、減振降噪、結構安全、尾跡隱身等密切相關，是船舶與海洋工程領域需要解決的重要問題。

當前，國內外針對船海結構物水氣泡混合流的研究熱點和難點問題主要包括2 個方面：

1） 水氣泡混合流的生成演化及其與船海結構物相互作用的機理問題。水氣泡混合流本身的物理機制高度復雜，涉及多相流、多尺度、劇烈湍流等問題的耦合，并且受多種影響因素的控制。考慮到其與船海結構物的相互作用，還需要探明水氣泡混合流與多系統運動邊界、繞流和分離流的耦合，這些問題將研究的復雜性提升到了更高水平。當前，機理探索方面的研究可以分為單純水氣泡混合流生成演化機理和船海結構物周圍水氣泡混合流機理兩方面。對于前者，需重點從微觀層面分析因多相流動界面撞擊作用而產生的氣泡和液滴結構，并從物理特性的角度研究水氣泡混合流的多因素影響規律；對于后者，則重點從宏觀層面分析船海結構物不同位置處水氣泡混合流的生成條件與演化分布，從流動、能量和載荷等方面研究水氣泡混合流與結構物性能的關系。機理問題的研究有助于深化人們對船海結構物水氣泡混合流問題的認識和理解，并為進一步對其進行控制和優化奠定理論基礎。

2） 水氣泡混合流的高可靠性數值模擬方法問題。水氣泡混合流的一個重要特點是其廣泛存在于實尺度結構物周圍，而在模型尺度試驗中則難以復現。在相似規律不明確的條件下，難以通過試驗手段對實際情況進行預報，制約了在船海結構物水氣泡混合流方面的研究突破。因此，發展高可靠性數值模擬方法就顯得尤為重要。當前，在數值方法方面的研究可以分為精細化模擬方法和模式化模擬方法。對于前者，重點是建立可高精度捕捉微尺度界面變形和氣泡、液滴結構生成演化的數值算法，解決界面捕捉、數值耗散、守恒格式等關鍵問題，用以進行船海結構物水氣泡混合流的微觀物理機制分析；對于后者，重點是建立能模擬實際大規模氣泡液滴群的精確物理模型，并與界面流動模擬方法相結合，形成跨尺度流動求解方法，用以對實際船海結構物周圍水氣泡混合流的宏觀特征進行預報。數值模擬方法的研究可以為流動機理分析提供有效手段，也能為開展實際結構物周圍水氣泡混合流動的預報和構型優化設計提供分析工具。

因此，本文將分別對機理探索研究和數值方法研究的進展予以介紹。

1 船海水氣泡混合流機理

1.1 水氣泡混合流生成演化機理

水氣泡混合流中包含著多尺度的流動結構[6]。在微觀上，有毫米及亞毫米尺度的微氣泡和微液滴；在宏觀上，有大氣泡、液體射流和大范圍的水氣自由面。微觀尺度的流動結構在表面張力的作用下可以保持一定的形狀，而宏觀尺度流動結構的特點是水氣界面存在連續不斷的變形。在生成演化機理探索上，目前國內外的研究可以分為4 個方面：一是以宏觀流動為基礎，研究局部區域內相界面的撞擊作用和破碎過程，并分析小尺度流動結構的來源；二是以微觀流動為基礎，研究氣泡、液滴結構的撞擊作用和演化形式；三是研究黏性和表面張力對水氣泡混合流的影響；四是研究水氣泡混合流與噪聲的相互作用。下面，將從這4 個方面分別介紹相關的研究進展。

1.1.1 局部區域多相界面撞擊作用

船海結構物周圍局部水體的水氣界面撞擊作用，是船海結構物水氣泡混合流生成的一個重要機理問題。在撞擊過程中，由界面破碎引起的非線性作用會產生氣泡摻混和液滴飛濺，過去，通常是利用波浪翻卷破碎試驗來模擬這種界面撞擊現象。在試驗中，宏觀的波浪破碎流動過程主要采用高速攝像系統進行記錄[9-11]。在早期研究中，通過現象觀察以及氣泡尺寸測量，發現水下氣泡群呈現明顯的多尺度特性[12]，如圖1 所示。氣泡直徑以Hinze 尺度為界，具有不同的生成機制，較大氣泡的生成機制屬湍流誘導破碎，其決定了水氣泡混合流中大尺度氣泡的體積分布；較小氣泡的生成機制主要是射流和液滴與自由面的撞擊作用。Kimmoun 和Branger[13]采用粒子圖像測速（particle image velocimetry，PIV）技術，通過測量水下速度場和漩渦流動，研究闡述了湍流是大氣泡破碎的主要原因，以及破碎后不同尺度的氣泡所導致的分布區別。直徑較小的氣泡受到的浮力較小，這些氣泡會在海洋中聚集與漂浮，形成氣泡層[9]。在自由面上方，波浪翻卷撞擊會在水面產生高速射流，并在空氣中分裂成大量的小水滴，然后在劇烈的流動條件下，水滴會在空氣中霧化，形成霧化水滴云。霧化水滴云和水下氣泡云是船海結構物水氣泡混合流的重要流動結構。

圖1 波浪破碎誘導的水氣泡混合流中的多尺度特性[12]Fig. 1 Multi-scale characteristics of water-air-bubble mixed flow induced by wave breaking[12]

由波浪翻卷撞擊破碎試驗，可獲得豐富的水氣泡混合流相關的流動現象，但試驗中最大的困難在于對微觀流動結構無干擾的定量統計和測量。為更加準確地探究波浪破碎的演化及破碎過程中氣泡、液滴運動的動力學特性，數值模擬方法被應用到波浪破碎問題研究中。早期，因計算條件所限，部分學者在湍流模擬上采用雷諾平均（RANS）或大渦模擬（LES）等模型化方法進行了模擬計算[14-16]，并結合了一些用于描述離散氣泡群效應的模型，重點是分析試驗中難以測量的流動速度場，分析波浪翻卷到波前產生飛濺的過程，并提取飛濺液滴以及氣泡摻混的情況，但該模型化方法的缺點是流場中的脈動并不能被準確求解，導致氣泡多尺度結構特征描述存在一定的誤差。

近年來，隨著計算技術的快速發展，許多學者[5,17-18]開展了波浪撞擊破碎過程的三維直接數值模擬（direct numerical simulation，DNS），典型的模擬結果如圖2 所示。從中可以看到，能精細捕捉到撞擊過程中的復雜界面破碎。借助DNS，還可以發現波浪破碎早期的渦絲現象，這也進一步解釋了波浪破碎過程中的卷氣機制以及相應的能量損耗機理。另一方面，DNS 的精細效果也推動研究者可以從宏觀?微觀的耦合層面來理解微氣泡等流動結構的生成。Chan 等[19]的研究表明，在一個27 cm 范圍內的水氣泡混合流中，每秒可發生上百次的微氣泡撞擊，可見，在實際中每立方米空間內每秒會產生上百萬個微氣泡。這種多尺度流動現象的探索對理解實際船海水氣泡混合流生成演化機理具有啟示作用。

圖2 波浪破碎中水氣泡混合流生成與渦絲的數值模擬結果[17]Fig. 2 Numerical simulation results of water-air-bubble mixed flow generation and vortex filament in wave breaking[17]

基于多尺度的精細化數值模擬，Wang 等[18,20]總結出了波浪破碎問題中的4 種水下氣泡產生機理，如圖3 所示，包括射流碰撞卷氣、碰撞形成二次飛濺卷氣、反向飛濺與初始波峰之間夾帶氣體（反向飛濺卷氣）和湍流夾帶卷氣。在不同的機理作用下會形成不同尺度的氣泡，而這些氣泡又會進一步呈現出不同的動力學行為。同時，水面以上霧化液滴的形成也被描述為液體表面拉伸、飛濺凸起和液體脫落。Wang 等[18,20]研究發現，微觀流動之間的相互作用發生在自由面附近，即氣泡在自由面附近破碎會生成新的液滴，而液滴下落砰擊自由面也會產生新的氣泡，這是除大尺度流動破碎演化之外微尺度流動結構新的來源。

圖3 破碎波水氣泡混合流中的4 種卷氣機制[18]Fig. 3 Four air entrainment mechanisms for water-air-bubble mixed flow of wave breaking[18]

1.1.2 微觀流動界面演化機理

船海結構物水氣泡混合流多相流動界面的相互撞擊過程實際包含了大量微觀流動結構的界面演化，如氣泡、液滴等。目前，一些研究人員主要從少量氣泡的相互作用，以及氣泡和液滴與自由液面相互作用的角度，對氣泡液滴生成微觀機理進行了研究分析，例如單個大氣泡的破碎過程試驗[21]、極少量氣泡間非線性相互作用試驗[10-11]、氣泡上升與自由面相互作用試驗[22]、液滴下落撞擊自由液面過程試驗[23]等。這些試驗研究利用高速攝相機記錄了水氣界面復雜演化的全過程，研究分析了氣泡破碎過程的發展規律以及末期射流的產生原因、氣泡間撞擊及聚并過程的非線性作用機理（圖4），為氣泡云及水氣泡混合流生成機制的研究奠定了基礎。

圖4 微觀界面撞擊形態關于Oh 數和We 數的相圖[23]Fig. 4 Regime maps for micro-scale interface morphology based on Oh and We[23]

在以上試驗中，測量設備和儀器會對氣泡、液滴的形成產生干擾。為克服該問題，近年來，研究人員開始采用計算流體力學（CFD）數值模擬方法開展研究。在氣泡的微觀界面撞擊作用機理方面，模擬研究了氣泡與氣泡之間的靠近追逐[24]、碰撞反彈[25]、接觸聚并[26]等過程。在靠近追逐問題上，共軸線氣泡上升是機理研究的經典問題，當前的研究已經涉及到了各種雷諾數范圍[27]，分析了上部氣泡的屏蔽作用，并在雷諾數Re≈ 100的情況下發現了穩定“平衡距離”現象[28]。然后以此為基礎，又進一步將研究拓展到了并排上升和多角度上升方面，其中尾流效應仍是研究的重點[29]。在接觸聚并問題上，目前有已得到廣泛認可的氣泡/液滴聚并過程描述[30]。對于氣泡來說，首先兩氣泡彼此靠近，在二者之間形成一層非常薄的液體層，然后隨著壓縮作用液體的排出，在達到臨界厚度時破碎，從而產生局部的微觀氣橋現象，最后兩氣泡中的氣體相互連通，即完成聚并。對于液滴，也有相似的氣體層和液橋現象[31]，主導著聚并過程的演化。聚并轉化的條件受表面張力，也即韋伯數We的控制。

在當前基于網格類方法的數值模擬中，大多是采用流體體積（volume of fluid，VOF）法或是水平集（Level-set）法進行水氣界面的捕捉，但這兩種方法均存在一定的不足。VOF 法采用連續流體模型計算界面表面張力存在較大誤差，而Level-set 法在保證質量守恒上存在問題。因此，有學者耦合這兩種方法的優勢，采用CLSVOF 方法進行了比較復雜的氣泡/液滴微觀多相流動界面撞擊流動模擬，包括單液滴撞擊傾斜液膜飛濺過程模擬[32]、三氣泡上升過程中的相互作用模擬[33]等。借助于該方法，可以更好地描述氣泡/液滴表面的變形和運動特性，從而能為進一步開展表面張力對流動特性的影響研究奠定基礎。另一方面，也有研究采用粒子法對高度分散的微尺度氣泡液滴生成演化過程進行了模擬[34-38]，并對翻卷破碎、射流等劇烈的自由面演化以及由其引起的流場脈動問題進行了分析。

1.1.3 黏性和表面張力影響機理

水氣泡混合流的生成演化受流體黏性、湍流流動以及界面之間表面張力等因素的影響，在有關影響因素的研究中，由于其參數調整的便捷性，理論研究和數值模擬成為主要的研究手段。同時，也有一部分試驗工作是通過采用大部分物理參數接近而某一物理參數不同的流體介質來對影響機理進行研究。

在有關黏性影響效果的研究中，Zhang 等[39-40]以不考慮黏性的勢流理論為基礎，通過引入黏性修正壓力和正應力的黏性分量對非零黏性剪應力進行補償，分析了有、無黏性導致的區別。該研究分析了單個氣泡在浮力作用下的變形、破碎以及射流過程，發現在考慮黏性的情況下，氣泡射流速度是隨流體黏性的增加而衰減的，氣泡射流能量降低，流場損失的動能和勢能逐漸轉化為黏性耗散能。同時，還研究了黏性對單個氣泡和雙氣泡脈動特性的影響規律，發現流體黏性會抑制氣泡脈動。以上研究的特點是計算效率比較高，可以開展系列研究并得到宏觀的規律。王昭太等[41]則采用基于黏流CFD 的模擬方法，探討了流體黏性對氣泡上升運動特性的影響，研究認為流體黏性主要影響氣泡內的速度分布，阻力系數隨黏度比的增大而減小。除了針對單個氣泡/液滴動力學行為的研究，有學者還通過數值模擬分析了流體黏性對水躍（hydraulic jump）等復雜界面流動中水氣泡混合流氣泡摻混特性的影響[42-43]，結果表明，隨著黏性作用的減小，摻混區內最大空隙率和氣泡數密度呈單調增加的趨勢。黏性影響和湍流效應密不可分，Thoroddsen 等[44]利用PIV 技術，開展了湍流流動中氣泡運動變形的力學行為研究，并得到了指定區域的剪切速率、湍流動能和動能耗散等信息，發現氣泡在湍流場中的破碎演化過程是由外部因素（如湍流強度、剪切速率）和內部因素（連續相的流體動力學因素和相界面相互作用）共同作用的結果。Zhang 等[45-46]通過數值模擬分析了湍流邊界層與氣泡群的相互作用，并給出了湍流作用下的氣泡遷移運動軌跡與尺寸分布特性。

對于表面張力影響研究，在單獨氣泡/液滴動力學行為方面，艾旭鵬和倪寶玉[40]研究了表面張力（We數）對單個球狀脈動氣泡的影響，結果顯示表面張力會縮短氣泡的振蕩周期但幅值不變， 同時還會使勢能增加。Ma 等[47]研究了表面張力對從邊界注入的單個氣泡生成過程的影響，結果顯示表面張力的增加會導致氣泡生成平均周期增加、氣泡體積增加、氣泡生成頻率降低。在氣泡群動力學行為方面，Gourich 等[48-49]以氣泡羽流問題為基礎，分析了流體表面張力對氣含率、氣泡形態和尺寸分布的影響，認為表面張力的變化會明顯影響氣泡群中的聚合現象，氣泡聚并速率將隨表面張力的降低而降低，水、氣兩相體系的傳質速率也會明顯降低，并且流場中的氣泡體積分布更小。Chaumat 等[50]進一步研究了表面張力對氣體破碎的影響，發現表面張力減小后，會使水氣泡混合流在低流動速度下更容易發生失穩，進入液體流場中的氣體破碎得更快，這同樣也是氣泡尺寸分布減小的重要原因。對以上文獻中介紹的影響因素、研究對象和影響效果予以了總結，結果如表1 所示。

表1 黏性和表面張力對水氣泡混合流的影響效果Table 1 Effect of viscosity and surface tension on water-airbubble mixed flow

1.1.4 對噪聲的作用機理

水氣泡混合流動內部包含數目龐大的氣泡和液滴微尺度流動結構，致使整個流場存在高度的脈動。這一特點使得水氣泡混合流會對水面船舶噪聲的產生與輻射產生重要影響。

首先，水氣泡混合流中氣泡和液滴的動力學行為會對流場脈動造成影響，產生新的聲源。Song 等[51-52]利用水聽器對氣泡形態變化過程中周圍聲壓的變化進行檢測，并在試驗中分析了氣泡大小、直徑、潰滅方向等因素對噪聲輻射的影響。同時，也有學者認為液滴產生的噪聲本質上也來源于氣泡[53]。當液滴撞擊自由面時，引起的摻混氣泡受到局部壓力的作用會產生振蕩，進而輻射噪聲。在聲理論中，考慮將氣泡構造為聲源的方法可以分為線性理論和非線性理論。當前，主要的方法是線性理論，即只考慮球狀收縮膨脹的單極子源[54]。但近年來的實驗研究證明，對于復雜的流動，非線性聲源對水下噪聲的影響同樣不可忽視。

其次，氣泡的存在會對噪聲的傳播產生影響，這主要體現在聲速變化和聲功率損失2 個方面[55-56]。當噪聲在流體中傳播穿過氣泡時，氣泡會影響聲速并耗散能量[57]。早期，Commander 等[58]針對線性聲波穿過水氣泡的問題進行研究，提出了一套預測方法。 在該線性波理論的基礎上，Louisnard[59]開始關注非線性成分對聲波傳遞的影響，并研究發現，在氣泡流中，非線性波相比線性聲波衰減得更明顯。隨后，Fuster 等[60]在此基礎上又進一步考慮了質量轉換、氣泡融合等動力學行為，發現可以更精確地預測水氣混合流環境下的噪聲。

1.2 結構物周圍水氣泡混合流機理

船舶和海洋平臺作為典型的船海工程結構物，在實際作業過程中，其周圍廣泛存在著水氣泡混合流現象，不過機理和作用效果有所不同。船舶周圍的水氣泡混合流主要由其航行運動引起，并具有明顯的位置區域差異性，即在船艏部、舯部和艉部不同位置處生成的水氣泡混合流現象與機理均不同。而海洋平臺周圍的水氣泡混合流則主要源于惡劣海況下的波浪砰擊和甲板上浪。下面，將從這幾方面來介紹相關的研究進展。

1.2.1 船艏水氣泡混合流現象機理

船艏興波是伴隨船舶航行運動而產生的，在實尺度下一般不是規則理想的波動，而是會出現船體艏部與流體的劇烈撞擊，出現波浪破碎和水氣滲混，進而形成水氣泡混合流。因此，船艏部的興波破碎成為研究的首要問題。經典的試驗研究主要針對的是不同弗勞德數Fr下的船艏興波形態[61-63]。在試驗中，Fr的范圍一般在0.28～0.45之間，可通過提高航速來誘導產生興波破碎，然后采用自由液面可視化技術觀測船艏興波破碎產生的過程、艏波破碎流動結構以及艏波在船體表面上的分離點，還可以利用PIV 技術測量船模周圍的速度場和渦量場，分析在興波起始位置處出現的負渦區和在船邊界波峰處出現的正渦區。試驗研究指出，船艏的幾何構型特征會使興波向兩側翻卷，在高速情況下，興波的劇烈翻卷撞擊是產生船艏部水氣泡混合流的主要原因。

前期的實驗研究針對多種船艏幾何構型和航速，給出了興波破碎生成的水氣泡混合流現象的結果，再運用數值模擬方法，可以進一步定量研究船艏興波破碎產生的水氣泡混合流的物理量。由于需要考慮復雜幾何曲面結構物的影響，數值模擬研究的難度要比前文所介紹的純波浪破碎模擬大得多，因此前期主要的數值模擬研究在湍流模擬方法上是以RANS 為基礎[64-65]，Fr數最高可達0.62。通過模擬，可以成功復現興波在船艏側面的翻卷和射流現象，但平均化方法對于破碎的細節卻難以捕捉。近年來，基于分離渦模擬（DES）的湍流模擬方法成為船艏水氣泡混合流模擬的主要手段[66-69]，Fr數最高進一步提高到了1.2。典型的數值模擬結果與實船海試結果的對比如圖5 所示，圖中A，B，C，D 分別為4 個典型的現象對比區域。從中可以看出，由數值模擬可以得到船艏射流分離、自由面大變形、興波多次翻卷以及水氣摻混現象。這些數值研究是利用軸向渦量、速度和自由面表面拓撲的截面分析來解釋興波破碎的演化的。首先，液體在船艏表面爬高并向外側翻卷形成一次破碎，在高速情況下，射流在與自由面接觸后會發生二次甚至是多次反射，形成多次破碎，并在射流與自由面接觸的位置留下明顯的痕跡。通過對渦量場的分析，發現射流末端以不連續的速度與自由面重新接觸后，在該位置會形成一個強反旋轉渦對，從而進一步形成了大橫流梯度和新的射流。從現象完整性的角度來看，數值模擬研究可以用于分析討論船艏水氣泡混合流生成前的過程和機理，但還不能對破碎后形成的大范圍水氣泡混合流進行準確復現。另一方面，當前的試驗和數值模擬強調在高Fr數下才能出現破碎現象，但實尺度船舶在以低速航行時仍會出現水氣泡混合流，這說明該現象并不能完全滿足弗汝德數相似，其深層的影響因素和流動機理仍需進一步研究。

圖5 船艏興波破碎水氣泡混合流數值模擬與實船海試對比[67]Fig. 5 Comparison between numerical simulation and sea trial of water-air-bubble mixed flow induced by bow wave breaking[67]

波浪與船艏的砰擊作用是船艏生成水氣泡混合流的另一個重要原因。在波浪的作用下，船艏抬升和下落時與水體撞擊會產生強烈的砰擊霧化和水下氣泡摻混，從而形成大范圍的水氣泡混合流。近年來，Bodaghkhani 等[70-72]利用先進的測量設備針對該現象開展了模型試驗研究。在試驗中，基于規則幾何平板與復雜幾何船舶模型，對波浪砰擊過程進行了研究，并采用氣泡圖像測速儀（BIV）和數字粒子圖像測速法（DPIV）這2 種圖像測速技術研究了破碎波在沖擊平板的過程中產生的霧化現象。值得注意的是，該研究結合了實船海試中的一些流動現象，通過與模型試驗相配合，對砰擊霧化現象建立了新的認識。針對船艏砰擊形成霧化液滴的問題，已形成了一些有效的理論模型[71]，并已針對霧化液滴尺寸、速度和空間分布進行了建模，可用于流動現象的快速預報。

1.2.2 船舯水氣泡混合流現象機理

船舶舯部附近是船舶水氣泡混合流演化的重點區域，在船艏形成的水氣泡混合流會在船舶舯部附近形成一定的分布。前面已提到，在模型尺度下由船模試驗產生水氣泡混合流需要很高的航速，且生成的水氣泡混合流也很難演化覆蓋整個船體區域，因此難以從模型尺度的角度研究該演化問題。Johansen 等[1]開展了很有價值的實船觀測研究，在研究中對實尺度下“雅典娜II”號船周圍的氣泡流場進行了測量，提出了一種考慮實際探頭尺寸的氣泡測量新方法，并通過實測對流場中的氣體體積分數、氣泡速度和氣泡尺寸分布予以了統計。經研究發現，由于船體周圍的加速度，在船舯處演化的氣泡平均速度要略大于船速，并在邊界層內部逐漸減小；氣泡速度隨著深度和側面距離的增加而增加，并且靠近自由面附近氣泡的平均尺寸更大；當在不同航速、不同位置逆浪航行時，摻氣量均隨航速的增大而增大。另外，研究還針對湍流邊界層氣泡卷吸以及尾部氣泡流動等現象進行討論，得到了氣泡速度關于弦長的分布，稱為“弦長分布”。

Ma 等[72]基于實船觀測，提出了Phenomenological Air Entrainment 模型，用于預報水面船舶周圍的氣泡，得到了船體周圍氣相分數分布規律，如圖6（圖中，εg為氣泡源項強度）所示。同時，還討論了在直航和定常回轉運動中，夾帶氣流的位置與產生氣泡流間的定性差異，分析了船舶附體周圍的氣泡聚集現象，并將仿真結果與海上實驗觀測結果進行了比較，所得結論可為后續通過數值模擬進行機理探索奠定了基礎。船舶舯部附近的水氣泡混合流中的氣泡摻混量極大，采用界面捕捉類方法模擬計算的成本很高，因此對于船舶周圍氣泡流的數值模擬研究，主要采用以雙流體方法（two fluid model）為基礎的模型化方法[73]，以便于分析船舯水氣泡混合流的宏觀分布和演化過程。

圖6 船舯水下摻混氣的數值模擬結果[72]Fig. 6 Numerical simulation results of air entrainment at the midship positions[72]

1.2.3 船艉水氣泡混合流現象機理

船舶艉部水氣泡混合流最重要的宏觀表現是船后大范圍的白色泡沫尾跡，且該處的水氣泡混合流同時還受船舶尾部湍流和螺旋槳擾動的影響。由于該問題對于船舶的多種水動力性能和聲光特性具有重要影響，近年來得到了重點關注。對于船尾白色泡沫尾跡的生成機制，目前的研究認為其主要有2 個方面來源：一是航行過程中由船尾附近自由面流場翻卷引起的氣體卷吸[74]；二是由螺旋槳高速旋轉所產生的空化氣體[75]。自由面卷吸的氣體將以氣腔的形式進入水下，同時還會受到海水壓力和船尾分離漩渦流動的影響，導致其表面張力因難以克服流場作用力而發生破碎，最終成為微尺度氣泡。空化的氣泡則源自水中固有氣核在壓力低于飽和蒸氣壓條件下時的汽化，同樣表現為微尺度氣泡。這兩方面的氣泡會在船尾形成具有一定長度的氣幕帶，也即白色的泡沫尾跡[76]。

在實測研究方面，重點關注了大尺度結構和宏觀參數。Peltzer 等[77]在早期對高速船尾部的水氣泡混合流進行了海上觀測，研究提取了船舶尾流中不同位置的多個橫向剖面進行流場分析，發現在所有尾跡剖面中均存在明顯受到表面張力作用的帶狀流動結構，該特征在船尾后方數十公里處仍然存在。通過聲吶技術，可以對尾跡水氣泡混合流進行成像[78-79]，Soloviev 等[78]利用三維聲吶技術，以內河航道中的船舶為基礎開展了實船試驗，分析了因自由面破碎、螺旋槳擾動所產生水氣泡混合流現象間的區別，進一步闡述了船體周圍氣體卷吸的位置和過程，并將實船尾跡的流動結構分析精細化到厘米尺度。

船舶尾部的水氣泡混合流現象與船艉型線有較強的相關性。許多試驗研究通過高速方尾模型試驗研究了尾流的特性。Shen 等[80]利用DPIV 方法對幾種不同幾何形狀船模的尾流水氣泡混合流進行了測量研究，重點關注了自由面湍流問題，并在研究中通過數據處理得到了船尾后方自由面附近流動結構和湍流的統計特征，包括不同截面位置的平均速度和脈動速度場。同時，還開展了基于Orr-Sommerfeld 方程的流體穩定性理論分析，發現在實驗中觀察到的彎曲特征是由平均尾跡的固有不穩定性所引起，導致流動出現了非線性特征。美國海軍研究部基于一個方尾的大尺度船模開展了不同航速下的高速拖航試驗[81-82]，其利用高速攝影、水下攝影、孔隙率探針、激光雷達等先進測量設備，從含氣率、氣泡流動、數密度變化、時空分布等角度研究了水氣泡混合流的特性，并分析了在不同船速下的區別和影響，典型的試驗圖像如圖7 所示。

圖7 高速方尾船模試驗中的水氣泡混合流[82]Fig. 7 Water-air-bubble mixed flow in the experiment of highspeed square-tailed ship[82]

在數值模擬研究方面，目前已有模型化方法的應用，主要以RANS 或LES 湍流模擬方法與雙流體兩相流模擬方法的結合為主[83-84]。Castro 等[84]重點討論了尾流不同位置處氣泡體積分布、湍流與氣泡流的宏觀相互作用和不同航速的影響規律等，結果顯示氣泡體積與尾流流向位置之間的關系符合指數函數，與船尾相近位置處的氣泡體積更小。另外，研究還將尾流中水氣泡混合流的生成演化區分為了幾個獨立的過程，包括由自由面附近由單渦引起的氣泡卷吸、復數渦對氣泡卷吸的疊加作用、氣泡群尺寸的分布及水下擴散過程，完成了物理建模并與試驗結果進行了對比驗證。Hendrickson 等[85-86]利用隱式大渦模擬（ILES）和浸入邊界法（IBM）對方尾尾流水氣泡混合流進行了精細化的模擬，討論了半寬吃水比對尾流演化形態的影響規律，通過氣泡識別算法提取了不同情況下水氣泡混合流的氣泡尺寸、含氣率和分布規律，并以水躍問題為參照研究了速度場、兩相流分布和湍流特征間的關系。

1.2.4 平臺水氣泡混合流現象機理

海洋平臺因長期作業于固定海域上，容易受到惡劣海況的影響，而波浪砰擊和甲板上浪是海洋平臺周圍水氣泡混合流生成的主要原因。目前，相關的機理研究主要針對氣泡摻混對結構物所受流體載荷的影響。在大波陡畸形波浪砰擊海洋工程結構物的過程中，有一部分氣體會被卷入水中，形成多尺度氣泡摻混的水氣泡混合流，這些氣泡會導致脈動增強，從而進一步產生周期性的載荷。圖8 所示為波浪砰擊直墻結構的載荷測量結果[87]，由圖中可以明顯看出水氣泡混合流的影響，而采用傳統理論模型（藍色曲線）則無法預報該現象。Hu 等[7]在實驗中制造出了不同程度的破碎波，用以研究不同狀態水氣泡混合流砰擊海洋平臺壁面時產生的載荷特性，包括輕微破碎摻混、帶射流翻卷、含大尺度氣腔和完全破碎摻混4 種情況。結果表明，含有大尺度氣腔的砰擊過程在直墻結構上產生的載荷最大；在輕微破碎摻混情況下，波浪爬升的高度最高。同時，因存在氣泡摻混的情況，給載荷帶來了明顯的脈動。該研究說明了不同尺度氣泡在水氣泡混合流砰擊過程中的作用與區別。Sun 等[88]針對摻混氣泡波浪砰擊固壁表面的問題建立了邊界元分析模型，其中重點考慮了翻卷中心的大氣腔影響，忽略了邊界破碎的小氣泡影響，同時分析了氣腔體積對壓力的影響。目前，大多數相關問題的研究還是以考慮大氣腔為主，而采用模型試驗和數值計算均很難復現出實際高度破碎的多尺度水氣泡混合流現象[70]，因此對于其影響機理還需要進一步開展方法和理論研究。

圖8 氣泡摻混引起的砰擊載荷變化[87]Fig. 8 Effect of air entrainment on impact load[87]

除了由于波浪砰擊產生的水氣泡混合流外，有研究人員還重點關注了月池結構處產生的水氣泡混合流[72]。在當前大型探測船和海洋平臺結構物上的大尺寸月池中，自由面的演化過程十分復雜，產生的水氣泡混合流會在水下擴散，進而干擾水下設備的使用。Yao 等[89]針對該問題提出了一種抑制氣泡群產生和下掃的阻尼裝置，并且通過實驗分析了水氣泡混合流的流動特性和裝置的抑制效果。

2 水氣泡混合流數值方法

2.1 水氣泡混合流精細化模擬方法

船海水氣泡混合流現象中包含著復雜的相界面變形和氣泡、液滴生成演化物理機制問題，同時，還涉及到多尺度兩相流和黏性流、湍流、分離流的耦合作用[18]。為了探明其機理機制，有必要開展精細化數值模擬研究。在網格類方法中，DNS 方法是一種運用三維非穩態的N-S 方程對湍流進行直接數值計算的方法，它不需要引入任何湍流模型，可與高精度的相界面捕捉方法相結合，在研究復雜的兩相流動的基本物理過程中起著十分重要的作用。同時，其還能實現對微觀流動機制的分析，并對氣泡、液滴、湍流渦等難以解析的流動結構進行精準捕捉。DNS 模擬方法對數值格式分辨能力和網格分辨率的要求很高，在算法構建上涉及求解精度、格式耗散、壁面處理、計算效率等多方面的因素，接下來，將從自適應直角網格、相界面求解方法、大密度比兩相流求解方法、表面張力模型和浸入邊界法5 個核心技術方面介紹相關的研究進展。另一方面，無網格粒子法可以通過自由面粒子來表征飛濺、射流、氣泡、液滴等強非線性兩相流結構，由于不存在對流項離散耗散等問題，也可以作為水氣泡混合流的精細化模擬方法，下面將予以介紹。

2.1.1 自適應加密直角網格

采用傳統均勻直角網格開展兩相流問題的DNS 數值模擬，在追蹤小尺度非定常流動結構時對網格量的需求明顯增加，存在計算效率低的問題，為此，發展了自適應網格加密（adaptive mesh refinement，AMR）方法。AMR 通過采用一系列矩形（多層）網格嵌套組合來控制方程在離散求解過程中網格拓撲結構的動態變化，然后根據流場的特點及計算需要對局部區域內網格進行加密，故采用AMR 方法能夠在保證在局部網格分辨率較高的情況下提高整體計算效率。當前，已得到廣泛應用的AMR 方法主要包含以下3 種：

1） 基于網格點的自適應加密策略，如圖9 所示。該技術是以網格點作為樹結構的單元，然后針對每個網格點進行劃分。Popinet[90]所開發的軟件包Gerris 就是典型的基于網格點的AMR 技術，適用于復雜幾何外形和界面流動的數值模擬。這款軟件包最開始是被應用于具有非規則邊界的區域，例如湖泊的淺水方程數值模擬[91]。2009 年，Popinet[92]在Gerris 軟件包中植入了VOF 方法和連續表面張力（CSF）模型用來模擬多相界面流動，其利用基于網格點的AMR 網格，得到了一些高精度的結果，如液滴飛濺[93]、波浪破碎問題[11]等。該方法的優點是可針對每個網格單元進行細分，非常靈活，加密效率較高，適用于有限體積方法（FVM）；缺點是必須從不同空間分辨率的相鄰網格中取得數據，并行效率低，負載平衡困難。

圖9 基于網格點的自適應加密策略Fig. 9 Adaptive grid refinement method based on grid points

2） 基于網格塊的自適應加密策略，如圖10所示。該技術是以網格塊作為樹結構的單元，初始計算時，由一個或者多個網格塊組成整個計算區域，每個網格塊在每個方向上具有相同的固定數目的網格，且每個網格塊可以單獨進行求解。當1 個網格塊被加密時，它會被4 個（二維）或8 個（三維）子網格塊所替代，每個子網格塊具有和父網格塊相同數目的網格且子網格塊網格的大小只有父網格塊的一半。天體物理的求解程序Flash[94]是一種典型的應用塊自適應網格的加密技術。許多研究人員利用基于塊結構的AMR 技術求解了多相界面流動問題，包括射流霧化[95]、移動接觸線問題[96]等。基于網格塊的 AMR 技術的優點是并行效率高、負載平衡容易實現，并且由于每個網格塊都是由幾何結構完全一致的均勻網格組成，故易于實施高階差分方法；相對的缺點是其加密操作需要以塊為基本單位，計算域整體的加密效率較低，網格量較大。

圖10 基于網格塊的自適應加密策略Fig. 10 Adaptive grid refinement method based on grid blocks

3） 基于網格片的自適應加密策略[97]，如圖11所示。該方法是以網格片作為樹結構的單元，既不要求父網格片與子網格片覆蓋的范圍相同，也不要求每個網格片包含的格點數目相同。該方法綜合了前2 種方法的優點，計算區域中的每一個網格都需要單獨被檢查和操作，加密效率高，但缺點是網格片覆蓋的范圍變化不定，算法復雜，且網格周圍用于插值模板點的數據結構復雜，通信頻繁不利于程序并行化，負載平衡較難。

圖11 基于網格片的自適應加密策略Fig. 11 Adaptive grid refinement method based on grid patches

2.1.2 相界面求解方法

構造追蹤相界面運動的守恒格式，實現高精度的界面重構，是精確模擬水氣泡混合流的重要問題之一。采用數值模擬方法實現界面間斷關系的先決條件是準確獲得相界面的實時位置和形狀信息。由于相界面的位置和形狀是隨時間變化的，因此需要針對界面演化建立模型方程。如果在流動過程中兩相流體界面處不發生相變、滲透及相互溶解等現象，可以認為相界面是一個物質面，即界面始終由同一批質點組成。近年來，追蹤界面演化的數值模擬方法得到了快速發展，具有代表性的界面數值方法有陣面追蹤法（front tracking）[98]、Level-set 法[99-100]、VOF 法[101]和相場法（phase field）[102]。依據相界面的描述方式，界面數值模擬方法可以劃分為界面追蹤方法（interface tracking）和界面捕捉方法（interface capturing）兩類。

界面追蹤方法是直接利用物質面假設，采用拉格朗日粒子來描述界面。若要通過當地流場速度推進粒子在流場中運動，可以根據粒子軌跡得到界面隨時間的演化情況。陣面追蹤法屬于典型的界面追蹤方法，Unverdi 等[103]將陣面捕捉方法與界面追蹤技術耦合得到了一種新的陣面追蹤方法，該方法的原理為：一套網格用于求解流體流動；一組自適應元素用來標記運動界面。這樣，粒子運動的速度就可以通過歐拉網格的流體速度插值得到。由于所使用的拉格朗日粒子之間的間距一般比歐拉網格的更小，故其界面分辨率較高，可以準確捕捉相界面，但在氣泡運動中不能捕捉到氣泡的大變形和破碎，以致不能得到廣泛的應用。隨著粒子法的不斷發展，Leung 等[104]使用均勻無網格粒子對界面進行追蹤，將粒子與最近的網格進行了關聯，結果顯示采用該方法可以很好地處理/控制界面拓撲變化。

界面捕捉方法是先定義一個標量場作為相界面的描述函數，然后再采用該標量場的特定等勢線來描述界面位置。由此，界面追蹤問題即轉化為歐拉網格上標量場隨時間的演化問題。典型的界面捕捉方法包括Level-set 法、VOF 法和相場法，它們分別選取符號距離函數、流體體積分數和相場函數來作為相界面函數。在標量場時間推進過程中，其局部分布可能會嚴重偏離理想的分布，從而導致界面法線方向及界面曲率等一些重要界面信息的計算產生誤差甚至是錯誤，因此需要對標量場進行人工修正。不同的界面捕捉方法采取不同的措施來實現：Level-set 法采用重新初始化（reinitialization）方法[105]，VOF 法采用界面重構（interface reconstruction）方法[106]，相場法采用在演化方程中添加包含四階導數的小量函數方法[107]。VOF 法是一種比較成熟的界面捕捉方法，被廣泛應用于流體力學研究與實際工程領域，VOF 分數在網格內的分布示意圖如圖12（a）所示。基于幾何重構的VOF 法具備較好的質量守恒特性，但要實現三維相界面的高階幾何重構比較復雜，目前只有少量針對二維問題的基于最小二乘法的界面重構方法及二次樣條函數重構方法[108]。THINC（tangent of hyperbola interface capturing）方法[109]是一種新型的代數VOF 方法，無需界面重構，其特征相函數為雙曲正切函數，具體如式（1）所示。由圖12（圖中，i為網格點序號，C為網格中的體積分數，χ 為相函數）中的對比可以明顯看出THINC方法與一般VOF 法在界面表達上的區別，其能有效抑制由界面不連續引起的數值偽振蕩，通過結合權重線界面重構方法[110]，能夠達到與幾何VOF法相當的精度，且算法簡單，易于擴展至三維問題模擬中。

圖12 網格內相函數分布一維示意圖[110]Fig. 12 One-dimensional schematic diagram of the characteristic function in grid[110]

式中：H為特征相函數；x為場中某點的位置矢量；t為時間； ?為水平集函數； β為表達插值函數跳躍陡度的經驗系數。

采用VOF 法，利用體積組分來計算界面曲率（用于表面張力的計算）比較困難，而采用Level-set方法求解關于符號距離函數的輸運方程時，因符號距離函數在界面處連續光滑分布，使得一些關鍵界面幾何信息（最小距離、法線方向、曲率等）的獲取更為準確，不過傳統的Level-set 方法較難保持流體質量守恒。近年來，CLSVOF 法[111]得到了廣泛應用，其很好地結合了Level-set 方法中可以較為準確地捕捉界面曲率以及VOF 方法可以保證質量守恒的優點。盡管這在一定程度上增加了計算成本，但基于該方法求得的相界面曲率光滑性較好，便于實施表面張力平衡模型。

2.1.3 大密度比兩相流求解方法

船海結構物水氣泡混合流是由空氣和水組成的兩相體系，二者的密度比值較大，因而給精細化數值模擬帶來了一定的困難，特別是對于水氣泡混合流中的射流、飛濺等復雜界面變形，或是因強剪切作用而形成的兩相界面流動現象，在數值迭代過程中，在氣相中有可能會出現較大的偽速度，從而導致界面發生不真實的變形，致使空氣速度持續增長，加劇兩相流數值求解的不穩定性，甚至是引起迭代發散[112]。Ghods 等[113]指出，動量方程和VOF 方程離散方法的不同會導致動量和體積組分的輸運存在相位差，由此引起的界面位置誤差會導致在界面附近產生錯誤的動量輸運，而氣相因密度較小，容易形成較大的偽速度從而破壞數值計算的穩定性。采用守恒形式離散求解N-S 方程可在一定程度上抑制偽動量的產生。

基于該問題，有學者提出了不同方法用來提高大密度比兩相流動模擬的穩定性。Rudman[114]提出了一種從VOF 標量方程中構造通量密度的方法，其使用更加致密的網格來更新體積組分。首先，針對網格上VOF 相函數的分布曲線進行積分，計算得到了控制體邊界的相分數值；然后，通過插值得到了邊界的密度通量，以用于動量的輸運。該方法保證了動量和體積組分以相同的離散方法推進，實現了大密度比二維氣泡動力學行為的模擬。Bussmann 等[115]借鑒Rudman[114]的思想，將類似的算法移植到非結構化同位網格中，避免使用更致密的網格構造界面通量，降低了計算量，且應用于二維潰壩流動、大密度比液滴定速傳播等問題的模擬也驗證了算法的高魯棒性。Vaudor 等[116]修改了Rudman[114]的原始算法，其基于單網格系統建立了一種通過半網格質量守恒計算質量（密度）通量并用其推進動量的方法，如圖13所示。圖中，i為網格點序號，Δx和Δy分別為網格在x和y方向的尺度，u為流體速度，ρ 為流體密度，下標f 和c 分別代表面心值與體心值。可在每一時間層面的求解中通過引入密度輸運方程來更新密度，采用相同的質量通量推進動量，并在由密度輸運方程求解得到的偽密度場獲得下一時間步的速度場。因在動量及質量（密度）推進過程中使用的質量通量由相同的數值格式構造，故在每一時間步中能夠實現動量及質量的一致推進，則可避免因自由面位置誤差而引起的界面附近偽動量或偽速度的產生，提高了數值求解的穩定性。這種動量?質量一致推進的思想近期已被一些學者成功應用于大密度比液滴高速砰擊濕潤壁面[117]及液體金屬涂層過程[118]的模擬中，代表了大密度兩相流模擬方法的新進展和突破。

圖13 動量?質量一致推進法中的速度和偽密度控制體Fig. 13 Control volume of velocity and pseudo-density in the consistent mass and momentum flux computation method

2.1.4 表面張力模型

船海結構物水氣泡混合流中包含有許多微觀流動現象，而表面張力對這些微觀流動現象起著主導作用，對其中流動不穩定現象的發生和發展具有重要影響。水氣界面在表面張力的作用下會收縮并斷裂，從而形成小氣泡或小液滴。CSF 模型[119]是計算流體力學中常用的表面張力模型，其原理是將表面張力轉換為體積力然后加入N-S 方程中進行考慮。表面張力源項的計算式為

式中： σ為表面張力系數； α為流體體積分數； κ為界面曲率，其計算式為

式中，ni為相界面單位法向量。

傳統的表面張力模型在模擬流體的表面張力時，流體表面受到的表面張力都是法向力，很難保證流體表面的光滑性。汪歡歡等[120]在CSF 模型的基礎上對流體表面施加了切向力，以使流體表面更加光滑，結果顯示模擬的流體表面張力效果更好。

相比CSF 模型，尖銳界面力模型（SSF）[121]能夠在網格分辨率相對較低的情況下正確模擬該表面的張力效應。將SSF 模型與CLSVOF 方法及虛擬流體法（ghost fluid method, GFM）[122]配合使用，并將表面張力施加在兩相界面上，對于毛細不穩定性其有著較高的解析精度，因此可以實現對表面張力占主導作用的小尺度界面流動問題的精細化模擬。可采用CLSVOF 方法計算得到的符號距離函數 ?直接計算界面曲率κ：

其中，符號距離函數 ?經VOF 函數場重構[123]獲得：

式中：Ng為周圍交界面網格的數目；αg為網格g中的體積分數值；Dgij為網格g中界面線段至網格（i，j）的法向距離； θgij為網格g中界面線段中心至網格（i，j）中心連接線與界面法向的夾角。

完成曲率計算后，將表面張力以壓力梯度算子的形式加入N-S 方程右端，并通過GFM 法引入插值模板點，然后通過界面間斷關系構造面心壓力梯度，如圖14 所示。圖中：x為網格空間坐標；x=xin為相界面位置；xˉ為構造的虛擬網格點坐標；p為壓力；p1(x)，p2(x) 分別為氣相和液相的壓力。采用該方法，可以提高壓力泊松方程及速度修正方程中壓力梯度項的離散精度。

圖14 虛擬流體法構造的界面壓力梯度Fig. 14 Construction of the interface pressure gradient using ghost fluid method

2.1.5 浸入邊界法

船海結構物水氣泡混合流要考慮壁面邊界的作用。在兩相流精細化模擬當中，要計算固?液?氣三相相互作用的力學特性并模擬交界面的形態比較復雜，且對于船體等復雜曲面外形邊界的描述也比較困難，尤其是在考慮航行船舶等運動物體時，隨著邊界的運動，每一個時間步網格都需要進行更新以適應邊界的變化。為了解決這一問題，Peskin[124]提出使用固定笛卡爾網格的浸入邊界法，其基本思想是將復雜結構的邊界條件模型化為簡 N-S 方程中的體積力源項，使用簡單的直角網格系統不僅可以回避貼體網格生成困難的問題，還能極大地簡化網格生成工作，易于處理動邊界問題。

按照對體積力源項的處理方式，可以將浸入邊界法分為2 大類，即連續力法（continuous forcing approach）和離散力法（discrete forcing approach）。連續力法最開始主要用于處理彈性邊界問題，此時，體積力源項的表達與物面的材料性質有關，由力學關系可以直接得到其表達式。典型的代表是由Peskin[124]提出的浸入邊界法。隨后，Goldstein等[125]和Saiki 等[126]對該浸入邊界法進行了改進，其利用反饋力來表示固體邊界，通過指定流體速度在流場某一期望點為0 來滿足無滑移邊界條件。離散力法的體積力源項由離散后的控制方程計算得到，依賴于控制方程的離散形式，一般無法獲得其解析表達式，主要用于處理固體界面問題。典型的離散力法有直接作用力法[127]、切割網格法（cut-cell）[128]以及虛擬網格法（ghost cell）[129]。

浸入邊界法提出至今，已從計算精度、時間步長限制、高雷諾數流動等幾個方面進行了研究改進。Liu 和Hu[130]提出了一種適用于復雜運動物體周圍流場求解的浸入邊界法，即運用Level-set方法追蹤移動物面：

式中：點P和點A如圖15 中所示，其中點P為流場中某點，點A為距離點P最近的物面點，PA為兩點間的方向矢量；n為點P附近物面三角形面單元的外法線矢量，圖中n1和n2為n的示意舉例，nA為A點處的法向矢量。在該研究中，改進了速度邊界插值方式，用以抑制壓力振蕩，并允許在數值求解過程中使用較大的時間步長。隨后，又在前期工作的基礎上引入了基于移動最小二乘法的高精度邊界插值算法[130]，通過結合物體六自由度運動求解技術及網格自適應加密技術，實現了不可壓縮流體?剛體耦合運動的精確預報及流場求解。Zhang 等[131]將浸入邊界法與幾何重構VOF 法（VOF/PLIC）相結合，模擬了二維水平圓柱的入水過程，隨后，又進一步將二階精度的尖銳界面浸入邊界法應用于運動固體的模擬[132]中，研究了圓柱入水過程中自由面的變化。Yang和Stern[133]在考慮物面潤濕性影響的基礎上，將一種簡單的接觸角邊界條件與尖銳界面浸入邊界法相結合，模擬二維圓柱入水過程，結果顯示，引入接觸角邊界條件后，圓柱物面后方非物理的閉合空腔消失了，自由面形態更符合真實情況。

圖15 用Level-set 方法判斷的網格與物面相對位置[134]Fig. 15 Determination of the relative position of grid points and object surface using Level-set method[134]

2.1.6 無網格粒子法

以上介紹的方法均為針對水氣泡混合流精細化模擬的網格類方法，雖然采用了諸多技術來克服數值耗散和邊界處理等問題，但受限于歐拉求解框架的限制，在復雜運動邊界和大變形破碎自由面問題上仍有網格處理上的困難。與之相比，無網格粒子法是基于拉格朗日框架構建的求解算法，該方法將連續流體視為帶有速度、壓力等物理特性的流體粒子，基于梯度算子和拉普拉斯算子離散控制方程，可以有效避免對流項網格離散求解的耗散問題。同時，無網格粒子法可以通過判斷自由表面的粒子來方便地描述自由面的演化，避免了通過網格來描述自由面時的困難[38]。因此，近年來有學者開始采用無網格粒子法開展水氣泡混合流相關問題的模擬探索。

文瀟[135]發展了一種高精度和高穩定性的多相流移動粒子半隱式（parallelization of moving particle semi-implicit, MPS）方法，其采用各相同時求解的耦合策略，通過引入密度光滑格式、相互作用黏度、表面張力模型以及兩相粒子碰撞模型，實現了對復雜多相流動的模擬。然后，在此基礎上又面向大密度比的水氣泡流動問題，加入了基于粒子間平均密度的壓力泊松方程求解、適用于高密度比的壓力梯度模型，以及可壓縮?不可壓縮耦合模型，使得采用該方法能夠對水氣高度摻混時的界面壓力和速度實現準確求解。Wen 等[34]和陳翔等[38]基于該無網格粒子法，對高度分散的微尺度氣泡液滴生成演化過程進行了數值模擬，在模擬過程中，采用流體粒子可以很好地模擬相界面翻卷破碎、水躍、射流飛濺、氣泡群擴散等強非線性現象，且研究還進一步對因多相摻混引起的流場速度和壓力脈動進行了分析。為了加速計算，無網格粒子法可與GPU 計算技術相結合[36]，以克服因流體粒子數目過多帶來的計算效率下降的問題。

采用無網格粒子法模擬船海結構物水氣泡混合流問題的局限性主要包括2 個方面[136]。第1 個方面是計算效率的問題。對于三維物理問題的模擬，鄰居粒子數目的增加會明顯增加相互作用搜索和模型處理的時間，同時也會增加壓力泊松方程求解的維度與困難程度，致使計算效率明顯降低，導致該方法在三維大尺度問題上的應用出現困難。第2 個方面是邊界層黏性流動的求解問題。無網格粒子法很難像網格類方法一樣布置出致密的邊界層網格，采用該方法模擬高雷諾數黏性流動效應和流動分離效應存在一定的缺陷，其對實際的復雜船舶和海洋平臺周圍流場的模擬存在局限性。

2.2 水氣泡混合流模式化模擬方法

船海結構物水氣泡混合流問題涉及超大時空尺度跨越，其中微尺度氣泡液滴為毫米量級，聚并破碎等動力學行為只需要幾毫秒，而宏觀尺度水氣泡混合流尾跡卻可達千米量級，完全生成或是消除需要幾十分鐘[137]。同時，水氣泡混合流廣泛存在于實尺度結構物周圍，而在模型尺度結構物實驗中卻不明顯，因而有必要開展大尺度船舶海試和模擬研究。基于以上兩方面的原因，要針對實際問題進行完全精細化的求解非常困難，因此有必要對其中的大規模微氣泡群等小尺度結構進行模式化建模，先通過氣泡卷吸模型等描述其生成特性，然后再通過宏觀界面捕捉方法與分散兩相流模擬方法的結合形成跨尺度耦合方法，以對船海結構物水氣泡混合流進行模擬。下面，將從3 個方面介紹其數值方法的研究進展。

2.2.1 氣泡卷吸物理模型

氣泡卷吸是氣泡混合流問題中的一個重要現象，尋求正確的模式化建模方法來模擬氣泡卷吸現象，對于研究上述問題十分關鍵，同時，也是進一步研究氣泡混合流與結構物相互作用的基礎。

誘導氣泡卷吸現象的物理過程主要發生在自由面附近，根據Castro 等[84]和Li[138]的分析，該過程可以分為直接卷吸和渦卷吸2 種。其中直接卷吸包括因波浪破碎而產生的射流沖擊水面卷吸氣泡入水、由波浪翻卷形成的氣穴卷吸氣泡入水以及由飛濺液滴撞擊水面捕獲氣泡導致的氣泡卷吸等情況，分別如圖16（a）～圖16（c）所示。對于渦卷吸，因湍流以及其他因素導致的旋渦會使得自由面發生變形，當變形過大時，可能會誘導氣泡卷吸現象，氣泡此時受到湍流旋渦作用力、浮力和表面張力的聯合作用。通常，切向旋渦和垂向旋渦也會誘導不同的氣泡卷吸運動，分別如圖16（d）和圖16（e）所示，相比切向旋渦，垂向旋渦通常卷吸氣泡的入水深度更深。圖中，ω 為漩渦旋轉的角速度矢量。

圖16 誘導產生氣泡卷吸的物理現象[138]Fig. 16 Physical phenomena that induce bubble entrainment[138]

近年來，學者們提出了許多氣泡卷吸模型。Baldy[139]基于相似性原理，對自由面處因湍流誘導波浪破碎產生的氣泡卷吸現象進行研究，提出了相應的氣泡卷吸模型。該模型基于Clift 等[140]的研究，假設氣泡卷吸速率與湍流耗散率ε呈線性關系，研究得到在自由面附近，氣泡的聚集程度正比氣泡直徑d的負二次方，該卷吸模型表達式為：

式中：S為卷吸量； εˉ為單位質量流體的湍流耗散率； ρ為密度； σ為表面張力系數；K1為模型系數。

Shi 等[15]對由波浪翻卷破碎誘導的氣泡卷吸進行了研究。其結合Deane 和 Stokes[12]在實驗中得到的氣泡尺寸分布譜，將自由面處的氣泡聚集程度與剪切作用建立了聯系。氣泡卷吸源項模型可以寫為：

式中： εn,i為在氣泡數密度對流擴散方程中添加的源項； αb為氣體體積分數；Pr=μt|S|2，為剪切生成項，其中 μt為渦黏系數；Di為與氣泡直徑相關的概率函數；drb,i為氣泡半徑空間分布。

Shi 等[15]將該氣泡卷吸模型應用到了雙流體方法當中，結果顯示，該模型在參數選取合適的條件下可以捕捉到氣泡卷吸過程的主要特點，但對于波浪剛開始破碎時，存在較大氣穴，模型未能表現出好的模擬效果。

Ma 等[141]基于Baldy[139]的研究，同樣假設氣泡卷吸速率與湍流耗散率ε呈線性關系，基于Carrica等[142]提出的多分散雙流體方法，建立了氣泡卷吸模型，該模型可以寫為：

式中：S(ri)為單位體積產生氣泡的速率；f(ri)為氣泡群粒徑譜；Cb為控制卷吸氣泡量的系數； ρl為液體密度； αl為液體體積分數；ri為第i組氣泡的半徑；NG為劃分的氣泡種群數目。

Ma 等[141]采用求解器TRUCHAS，對波浪沖擊海灘翻卷破碎時的氣泡卷吸現象進行了模擬，發現氣泡尺寸的分布規律與Deane 和 Stokes[12]的實驗結果吻合較好，研究顯示，直徑小于1 mm 的氣泡是由射流與自由面的相互作用產生的，而直徑大于1 mm 的氣泡則主要是由氣穴破碎所導致，并且氣泡的存在會在一定程度上抑制由波浪破碎導致的湍流。

式（9）所示的氣泡卷吸模型在Derakhti 和Kirby[16]的研究中得到了進一步應用。Derakhti 和Kirby[16]在Carrica 等[142]和Ma 等[141]研究的基礎上，將VOF 方法和雙流體方法相結合，針對深水孤立波波浪破碎過程開展了數值模擬。其采用VOF方法捕捉自由面，雙流體方法模擬分散氣泡，湍流模型則使用 LES 方法，并考慮氣泡對能量耗散和水相與空氣相間動量輸運的影響，模擬得到波浪破碎時氣泡摻混結果與試驗結果[143]的對比如圖17 所示。圖中，x*，z*分別為根據波長無量綱化后在波浪傳播方向和水深方向上的距離，自上而下各圖分別表示在參考時間（t-tb）為 0.05，0.15，0.30，0.55，0.75，0.95，1.15，1.45 時的結果。從圖中可以看到，水氣泡混合流的時空分布得到了較好的復現。研究發現，在波浪破碎區域，氣泡誘導的能量耗散占總能量耗散的50%以上。

圖17 深水孤立波破碎數值模擬與試驗對比[143]Fig. 17 Comparison between numerical results and experimental results of deep-water solitary wave breaking[143]

Ma 等[144]進一步修改了該模型，其假設氣泡卷吸速率與流體局部湍動能和流體速度梯度的乘積成正比，提出的亞格子氣泡卷吸模型如下：

式中：S(x)為單位時間單位體積內引入的氣泡總數密度；g為 由重力引起的加速度值；cent為該模型系數； ?ent為界面厚度；un為流體速度在交界面法向處的分量；k(x)為某位置處的局部湍動能。Ma 等[144]采用雙流體方法對俯沖射流和水躍這2 個涉及自由面的氣泡卷吸現象進行模擬，得到了較好的模擬結果。

目前，大多數氣泡卷吸模型都是基于以前的經驗和實驗結果，對初始時刻和物理邊界處的含氣率與氣泡尺寸分布給予了一定的假設，并將氣泡卷吸速率和湍流進行了聯系以進行進一步的分析研究。對于不存在結構物壁面而只考慮自由面的氣泡卷吸問題時，例如波浪破碎，這些模型體現出了較好的模擬效果；但當考慮壁面存在時，由于結構物壁面誘導的湍流往往遠離自由面，其誘導的氣泡卷吸問題有所不同，因而以上模型將不適用[84]。

船舶結構物在海上航行時會卷吸氣泡，進而形成水氣泡混合流。由于有關船海結構物的氣泡卷吸問題涉及的尺度較大，采用直接模擬的方法需要較多的計算資源，因此有必要尋求合適的氣泡卷吸模型。目前，多是尋求能較好地反映氣泡卷吸物理特性的亞格子模型來對氣泡卷吸問題進行數值模擬。Moraga 等[145]根據流體局部速度是否超過閾值以及距離交界面的距離的準則，建立了可以確定氣泡卷吸發生位置的亞格子氣泡卷吸模型，其表達式為：

Moraga 等[145]采用多分散雙流體方法求解氣泡輸運，在一定程度上較好地模擬了DTMB 5 415船因艏波浪破碎導致的氣泡卷吸和“雅典娜”號船近船體尾跡中的氣泡卷吸現象。Ma 等[72]將文獻[144]建立的氣泡卷吸模型用于模擬“雅典娜”號直航和定常回轉時周圍的氣泡卷吸現象，并將結果與實驗結果進行了對比，結果表明采用該模型可較好地模擬船舶周圍流場中空氣體積分數的分布情況。Castro 等[84]針對旋渦與自由面的相互作用過程進行建模，建立了氣泡卷吸模型，并應用REX 流體力學求解器，采用多分散雙流體方法研究了實尺度下“雅典娜”號航行時船體周圍的氣泡卷吸現象，自由面采用單相Level-set 方法進行確定，其模型框架形式為：

該模型克服了之前模型中氣泡卷吸僅由湍流參數決定而與氣泡相對自由面距離無關的缺點，并充分考慮了氣泡卷吸過程，包括由單渦引起的氣泡卷吸、復數渦疊加作用、卷吸氣泡群尺寸和水下擴散演化過程等，這也進一步反映了氣泡卷吸的力學機理。通過與Tavakolinejad[146]的實驗結果進行對比，發現該模型可以較好地捕捉氣泡卷吸位置和卷吸程度。Li 等[147]基于Castro 等[84]建立的氣泡卷吸模型，采用多分散雙流體方法，結合重疊網格技術和Level-set 方法，研究了“雅典娜”號船在轉舵過程中流場和旋渦對氣泡卷吸的影響，以及氣泡輸運和尺寸分布特性等問題，得到了較好的模擬效果。

2.2.2 歐拉?拉格朗日跨尺度方法

采用歐拉?拉格朗日跨尺度耦合方法模擬水氣泡混合流的基本思想是：在歐拉框架下，采用界面捕捉方法模擬結構物周圍興波等自由面演化現象，在拉格朗日框架下，采用離子追蹤方法模擬微尺度氣泡的運動，忽略氣泡表面的變形。氣泡運動求解的控制方程為動力學方程，一般具有如式（13）所示的形式，其中待求解物理量為氣泡運動速度，方程右端為各種受力模型。這種耦合方法由于氣泡的描述不依賴網格分辨率，因此可以用于模擬亞網格尺度的氣泡群演化，相比精細化求解，計算效率明顯提高，并且離散相的求解不存在數值耗散，適于計算水下大規模氣泡群的體積、運動和分布。

式中：m為氣泡質量；v為氣泡運動速度；ω 為氣泡旋轉角速度；u為液體速度；ρl，ρg分別為液體密度和氣體的密度；FD，FL，FF，FB，FAM分別為阻力、升力、流體加速度力、浮力和虛擬質量力；CD，CL，CVM分別為阻力系數、升力系數和虛擬質量力系數。

目前，針對該耦合方法的研究大多處于算法的設計和優化階段，主要應用在一些基本流動的科學問題上，例如射流霧化、氣泡羽流等。Herrmann[148]將Level-set 方法與拉格朗日追蹤方法相結合構建了耦合方法，其中Level-set 方法可以用于模擬占據足夠多網格的連續相界面，通過基于尺度判別的轉化算法，流場中不能通過網格合理求解的小尺度兩相結構將自動轉化為拉格朗日粒子進行模擬。利用該方法，明顯提高了射流霧化過程的模擬效果[149]。Patkar 等[150]的耦合方法結合了獨特的轉化算法與拉格朗日撒點模型，并利用該耦合方法對許多科學問題進行了模擬，包括上升氣泡破碎、氣泡羽流、空化等，典型的計算結果如圖18所示。其中，（g）圖通過染色對流場中采用不同方法描述的氣泡予以了區分，從中可以看出，采用耦合方法可以有效解決因網格分辨率不足而引起的小尺度氣泡不能捕捉的問題。在歐拉?拉格朗日跨尺度耦合方法中，摻混的小尺度氣泡、液滴量一般不需要額外進行人工模型干預，但在Hsiao 等[151]的研究中，其通過精細的試驗測量和理論分析，建立了亞格子尺度摻混模型，并將其應用到了水平射流問題的耦合方法模擬中，結果顯示模擬結果與試驗結果吻合較好。

圖18 采用Level-set 和拉格朗日追蹤耦合方法模擬的多尺度氣泡流[150]Fig. 18 Simulation of multi-scale bubbly flow using Level-set/Lagrange coupled method[150]

空化是一類典型的多尺度兩相流問題，該現象常出現在船舶螺旋槳和舵等推進系統的作業過程中。Hsiao 等[149-152]在拉格朗日粒子追蹤方法中添加了基于Rayleigh-Plesset 方程的空泡求解算法，該方程將氣泡半徑作為未知量進行求解，可以用于模擬空化過程中氣泡的膨脹和潰滅，控制方程如式（14）所示。

式中：R，R0分別為t時刻和初始時刻的氣泡的半徑；us=uenc?ub, 其中ub為氣泡的移動速度，uenc為氣泡運動過程中周圍環境流體的運動速度；penc為周圍環境流體的壓力；pg0為初始時刻氣泡內的壓強；pv為蒸汽壓； γ為表面張力系數；μ為流體黏度系數；kc為多變壓縮常數。

Hsiao 等[149]基于該方法模擬了水翼和螺旋槳表面的空泡初生，尤其是葉梢部位，典型的模擬效果如圖19（圖中，P為水翼表面壓力）所示。由圖可見，水翼表面上的空泡初生為微尺度的小氣泡，該現象可采用拉格朗日點來表示；隨后，在發展的過程中發生膨脹和聚合，形成大尺度片空泡，該現象可采用界面捕捉方法表示；最后，脫落后再次演化形成擴散的小氣泡由拉格朗日點表示。結果顯示，利用耦合方法后整個空化過程的模擬效果得到了明顯提升。

圖19 采用Level-set 和拉格朗日追蹤耦合方法模擬的水翼空化流[149]Fig. 19 Simulation of hydrofoil cavitation flow using Level-set/Lagrange coupled method[149]

歐拉?拉格朗日跨尺度耦合方法的關鍵問題在于能否全面考慮多尺度流動現象之間的轉化。早期的一些研究提出的算法是半耦合的[153-154]，即只能考慮拉格朗日粒子在與自由面發生接觸時的聚并融合，而由自由面翻卷自然破碎形成的氣泡摻混卻未得到考慮。這種方法適用于一些流場輸入條件為微氣泡的流動，如氣泡羽流等，對于船海水氣泡混合流這種因界面自然演化而形成的多尺度流動則難以模擬。在后續的研究中，又建立了全耦合的歐拉?拉格朗日跨尺度耦合方法。Tomar等[155]通過采用識別計算域中空相分數不為0 的網格，判斷了界面模擬至拉格朗日粒子的轉化，同時，拉格朗日粒子至界面模擬的轉化可以通過二者間的間距來判斷。該方法思路簡單、直接，為后續研究奠定了基礎。隨后，又有研究人員在此基礎上從轉化判據的角度對多尺度耦合方法進行了優化。在Ling 等[156]構建的數值方法中，包含超過43個網格，且與自由面距離超過1 倍直徑的氣泡或液滴將被轉化為拉格朗日粒子。反之，與自由面距離小于1 倍直徑的拉格朗日粒子將被轉化為界面捕捉方法來求解。在Zuzio 等[157]的研究中，轉化氣泡所需包含的臨界網格數為83個，其采用自適應網格技術，通過結合拉格朗日模式化方法，有效減小了自適應網格的加密等級，擴大了時間步長，提高了計算效率。需要注意的是，以上研究中轉化算法的實施需要以氣泡尺寸識別為基礎，常用的識別算法為泛洪填充算法，涉及到對全場網格的遍歷和臨近單元的搜索、識別與儲存。由于此過程在每個時間步都要執行一次，因而會增加計算時間并給內存帶來壓力，一定程度上限制了計算效率的進一步提升。Zhang 等[158]針對該問題提出了一種基于曲率的轉化算法，通過標模計算，對比說明了該方法在計算精度和效率上的優勢。

目前，歐拉?拉格朗日跨尺度耦合方法因兼顧了流動捕捉效果和計算效率的優點，所以在一些小規模的多尺度兩相流問題中得到了應用，但在船海工程領域應用得較少。對于實際的船海工程水氣泡混合流，微尺度氣泡液滴結構的數目在百萬、千萬級以上，而采用拉格朗日追蹤方法，由于存在粒子相互作用碰撞搜索、并行效率等算法方面的問題，導致在離散相數目規模擴大后計算量有明顯的增加，也很難用于處理大規模的船海工程水氣泡混合流問題，因此仍需進一步采用模式化方法。

2.2.3 歐拉?歐拉跨尺度方法

采用歐拉?歐拉跨尺度耦合方法模擬水氣泡混合流問題的基本思想是：對于大尺度的自由面，采用VOF 或Level-set 等界面捕捉方法進行捕捉，而對水下摻混的氣泡等小尺度界面，則采用雙流體模型進行求解，其中雙流體方法在求解過程中將氣泡混合流分為了連續相（液體）和離散相（氣泡），并采用歐拉方法進行了建模。求解的有關連續相和離散相各自的質量守恒及動量守恒方程為[145-146]：

式中：下標k代表連續相和離散相，取不同的值時代表不同的相；Mk為連續相和離散相之間的動量輸運，包括虛擬質量力項、拖曳力項、升力項、壓力項和浮力項等；Tk為應力張量， αk為體積分數；ρk為密度；uk為速度；g為重力加速度。歐拉?歐拉跨尺度耦合方法對于大尺度界面和小尺度界面的處理均采用了歐拉觀點。

目前，對于該種方法的耦合主要存在2 個方面的難點[159]。首先，界面捕捉方法與雙流體方法的求解框架不同。對于界面捕捉方法，兩相是通過求解一組動量方程來獲得速度場，因此速度場是共享的，但對于雙流體方法，速度場是通過求解兩相各自的動量方程組得到的，速度場是不共享的。其次，對相分數的理解這2 種方法也不同，對于界面捕捉方法，為了捕捉較為精確的界面，體積分數的擴散效應很小，但是對于雙流體方法，體積分數的求解強調了網格的平均效應，擴散效應會更加明顯。

歐拉?歐拉跨尺度耦合方法在波浪破碎、液體射流和水躍等不含結構物的水氣泡混合流問題的研究中已得到應用。Shi 等[15]基于Buscaglia 等[160]研究的雙流體模型和Deane 和Stokes[12]的實驗數據，發展了一種VOF 方法與多分散雙流體方法相結合的耦合方法，即采用VOF 方法模擬氣液界面，采用基于求解氣泡數密度的多分散雙流體方法模擬氣泡群，來對波浪破碎問題進行研究。Yan 和Che[159]建立了關于VOF 方法和雙流體方法的耦合模型，如圖20 所示，其將所有流體分為了連續液體、具有大尺度界面的大氣泡和具有小尺度界面的離散氣泡三相，這三相的體積分數分別為 α1， α2， α3，提出采用“體積分數再分配”的方法處理純流體、大尺度界面以及小尺度界面在網格中共存時產生的矛盾。經研究發現，對于小尺度流動的捕捉，該方法的求解效果與雙流體方法一致，而對大尺度界面的捕捉，該方法與VOF 方法的求解效果一致。

圖20 多尺度方法中流體相的劃分[159]Fig. 20 Division of fluid phases in multi-scale method[159]

對于垂直的液體射流撞擊自由面和水平水躍的問題，是不考慮結構物的典型水氣泡混合流的物理問題的。Xiang 等[161]建立了VOF 方法與雙流體方法的耦合模型，其中雙流體方法通過結合多種群方法求解氣泡尺寸，研究了水躍問題過程中的氣泡卷吸以及氣泡尺寸分布的變化情況，其與實驗數據[162]的對比顯示吻合較好。Ma 等[144]結合雙流體方法與Level-set 界面捕捉方法，提出了氣泡卷吸模型，并針對垂直液體射流撞擊自由面和水躍的問題予以了研究，如圖21 所示。圖中，α 代表空氣的體積分數，X為水流流向方向，Y為水深方向，X和Y均通過上游水深Dj進行了無量綱化。模擬結果與實驗結果的對比顯示，建立的模型可以較好地模擬相分數分布。

圖21 歐拉?歐拉跨尺度耦合方法模擬不含結構物的水氣泡混合流問題[144]Fig. 21 Simulation of water-air-bubble mixed flow without structure using Euler-Euler multi-scale coupling method[144]

在歐拉?歐拉跨尺度方法當中，對摻混氣泡群的另一種模式化求解思路是以氣泡數密度為基礎構建控制方程。基于玻爾茲曼理論得到的離散相氣泡數密度方程形式如下[142,163-164]：

式中：下標g代 表不同的氣泡種群；ug為氣泡速度；Ng為氣泡數密度，可用于進一步積分獲得氣泡體積分數； υt為無量綱渦黏系數；S cb為氣泡的施密特數；Dp為僅與體積分數有關的擴散系數； βg， χg，Sg分別為氣泡因破碎、聚并和卷吸而導致的生成源項。

由下式，可以獲得氣泡質量的輸運，同時，還需要配合描述氣泡速度的方程來表達動量守恒。動量方程的形式如下[142,163-164]：

在此基礎上結合基于Rayleigh-Plesset 方程的質量輸運模型，還可以實現空化過程的模擬[164]。

基于該思想，Carrica 等[142]首先建立了用于研究船舶水氣泡流動問題的歐拉?歐拉跨尺度耦合方法求解框架，其采用多分散雙流體方法并結合Level-set 方法，對波浪中水面船舶周圍不可壓縮水氣泡流動問題予以了研究。Castro 等[84]在Carria等[142]研究的基礎上，發展引入了氣泡卷吸模型，考慮到螺旋槳、舵以及其他附體等的存在，結合動態重疊網格技術，研究了“雅典娜”號船在直航情況所產生的氣泡混合流問題。Li 等[163]提出了一種可以改善求解氣泡混合流問題的魯棒性方法，改善了壓力速度耦合算法，可以實現高含氣率下的水氣泡混合流問題，成功對實尺度下一艘平底船周圍的水氣泡混合流現象予以了模擬。Li 和Carrica[165]考慮到氣泡的可壓縮性，提出了可壓縮氣泡混合流問題的求解方法，并模擬了“雅典娜”號船在2 個不同航速下的氣泡卷吸現象，結果發現船舶航行中，在船體表面、船艉以及船艏波浪破碎處會有較多的氣泡混合流現象發生。相對于高速情況，在低速情況下艉部的氣泡卷吸將更加明顯，含氣率也會更高。Li 等[147]基于以上研究方法，將其應用到了實尺度“雅典娜”號船在Z 形操縱運動下的水氣泡卷吸現象研究中，如圖22（圖中，Sg為氣泡卷吸源項大小）所示，得到了其在操縱過程中船體周圍水氣泡混合流的密度、含氣率、渦結構等特征，研究發現在操縱過程中，船艉處的氣泡卷吸程度會加劇，但在船體表面會減弱，其尾跡相對直航情況下也會更寬。

圖22 歐拉?歐拉跨尺度耦合方法模擬“雅典娜”號船在操縱運動下的艏艉含氣率分布[147]Fig. 22 Simulation of bow and stern void fraction distribution of Athena ship in maneuvering using Euler-Euler multi-scale coupling method[147]

相對于直接數值模擬方法和歐拉?拉格朗日跨尺度耦合方法，由于引入了更加模式化的方法求解氣泡群，歐拉?歐拉跨尺度耦合方法屬于多相流問題中的宏觀建模方法，在求解過程中不對氣泡的具體界面和形狀進行捕捉，而是著眼于捕捉流場中離散相的相分數分布。因此，此方法的求解效率會更高，對于求解實際船海結構物水氣泡混合流這種跨尺度問題的研究具有潛在優勢。

3 總結與展望

水氣泡混合流是一類在船海結構物周圍廣泛存在的復雜流動現象，其包含水、氣、泡、霧等多種流動狀態，具有多相、多尺度、多物理場等特點。目前，主要的研究重點在流動機理和數值模擬方法兩方面，本文以這2 個核心問題為主線，總結凝練了當前關鍵的研究方向和研究進展，目標是為研究人員認識、理解船海結構物水氣泡混合流機理提供基礎，為進一步形成有效的預報、分析和設計技術提供參考。

在流動機理方面，以從微觀到宏觀的角度系統介紹了目前船海水氣泡混合流機理的研究現狀。在微觀上，總結了當前以波浪破碎和氣泡/液滴等微觀界面流動為對象的水氣泡混合流生成演化機理探索工作，從局部水體多相界面撞擊角度介紹了自由面破碎和氣泡、液滴結構的生成機制，進一步介紹了生成后的微觀界面演化機理研究。從黏性和表面張力這2 個方面闡述了當前對水氣泡混合流動影響因素和影響效果的討論。最后，介紹了水氣泡混合流對噪聲源和聲波輻射的影響研究。在宏觀上，以船體和平臺周圍不同位置流動差異為線索，分類闡述了水氣泡混合流現象機理及其對結構物性能影響的現有結論，包括船艏、船舯、船艉型線導致的水氣摻混機理差異和演化區別，介紹了水氣泡混合流對平臺砰擊和上浪性能的影響研究。這有助于研究人員更好地認識和理解該現象，拓展研究思路。

在數值模擬方法方面，從精細化和模式化2 個角度，詳細介紹了適用于水氣泡混合流研究的數值模擬方法。在精細化方法上，首先針對網格類方法從求解精度、效率、穩定性等幾個方面介紹了AMR 方法、質量/動量一致推進、尖銳界面力模型、浸入邊界法等技術手段。對于無網格粒子法，介紹了多相流MPS 方法在劇烈自由表面流動和水氣泡混合流上的研究應用。在模式化方法上，介紹了當前用于模式化大規模微觀結構的物理模型以及基于歐拉?歐拉觀點和歐拉?拉格朗日觀點的模式化方法。在以上數值方法介紹中，分析了當前算法的優、缺點和能力，為研究人員進一步開展數值方法研究奠定了基礎。

隨著船海結構物不斷朝著高速化、大型化的趨勢發展，水氣泡混合流的影響問題也日益凸顯，在綠色節能、穩定安全、降噪隱身等方面體現出越來越強的研究需求，對未來的研究展望可總結為以下幾點：

1） 針對實尺度船海結構物的水氣泡混合流高效、可靠模型化預報方法。根據現有的觀測研究，船海水氣泡混合流的一個重要特點是在實尺度結構物周圍廣泛存在而在模型尺度卻難以出現。在相似規律不明確的現狀下，亟需發展針對實尺度船海結構物的水氣泡混合流高效、可靠模型化預報方法，建立能反映真實流動效應的物理模型，解決計算量、預報精度等關鍵問題，為真實流動分析提供手段。

2） 水氣泡混合流對船舶快速性能的影響機理。目前，船海結構物水氣泡混合流機理研究的重點在于對流動形態和生成演化過程的分析，對船舶具體性能影響的研究較少。當前，已有部分研究針對水氣泡混合流生成演化過程中的能量傳遞進行分析，但其消耗船舶能量以及引起阻力變化的定量化結論還有待闡明。另一方面，船舶尾流和船側的氣泡下掃可能導致氣泡群進入伴流區域，需要系統研究氣泡下掃量以及其對螺旋槳推進、空化性能的影響機理。

3） 水氣泡混合流對結構物砰擊上浪性能的影響機理。結構物砰擊過程本身由于瞬時性、極端性等特點具有一定的研究難度，而對純水波浪砰擊和含有大氣腔的波浪砰擊問題則有一定的研究。根據實際觀測，水氣泡混合流是結構物砰擊上浪過程的重要流動形式，需要進一步探明多尺度氣泡、液滴群對砰擊載荷上浪載荷所產生的非線性影響，研究其引起的砰擊上浪范圍、幅值的變化特性。

4） 水氣泡混合流與噪聲的多物理場耦合機理。氣泡、液滴的振蕩、聚并和破碎行為以及其與自由面的撞擊等可被視為重要的噪聲源。船海結構物水氣泡混合流中包含數目龐大的氣泡液滴結構，對水面艦船流動噪聲的影響是未來的一個重要研究方向，需要構建兩相噪聲預報技術，研究水氣泡混合流產生的船舶近場噪聲源分布，并在此基礎上分析遠場聲輻射特性，掌握水氣泡混合流與噪聲的耦合機理。