駐留式微氣泡陣列流動減阻機理數值研究

朱睿，何星宇，趙晨鴻，劉宇，張煥彬，陳騰飛，譚鑫，劉志榮

(1.廈門大學 航空航天學院，福建 廈門 361005；2.西藏民族大學 信息工程學院，陜西 咸陽 712082)

水下航行體的阻力主要有興波阻力、壓差阻力和摩擦阻力[1].調整形體和流線型可以改善興波阻力和壓差阻力，而低速時航行體的摩擦阻力約占總阻力的80%～90%[2].隨著數值計算和科學試驗的發展，對摩擦阻力的機理研究愈發深入.改變流體介質運動特性和邊界層內近壁區流場特性可以大幅減小航行體的摩擦阻力，是未來流動減阻技術的發展方向.

微氣泡減阻技術的試驗和數值結果顯示其良好的性能和應用潛力，已成為水下減阻領域的研究熱點[3-7].Kim 等[8]研究水翼表面在不同弗勞德數下的噴氣減阻特性，利用高速相機量化了氣泡注入量，取得了一定的減阻效果.Mohammed 等[9]采用數值方法研究集裝箱船模型在不同工況下的噴氣微氣泡減阻效果，結果表明當弗勞德數為0.282，空氣體積分數為4.8%時，最高減阻率達到27.6%.Jha 等[10]測試了不同氣泡空隙率α=0～0.15、雷諾數Re=22 500～67 500 及氣泡噴射方向下水平紊流單通道的流動阻力，結果顯示，頂壁噴射能夠減少最高60%的阻力.陳正云等[11]通過試驗觀測，得到不同超疏水表面形成的表面氣膜具備良好的穩定性和減阻性能，減阻率可達30%.姚琰等[12]通過試驗測量了多孔噴氣平板表面在不同通氣量和來流速度下的切應力，結果顯示微氣泡在浮力作用下形成邊界層氣膜，最大減阻率可達36%.宋武超等[13]通過試驗研究不同俯仰狀態下航行體表面的微氣泡形態及減阻特性，結果顯示，微氣泡附著受通氣量、重力和來流速度的影響，減阻率隨著航行體攻角的變化呈現正弦變化.

當前，主流微氣泡減阻方式通常是通過近壁區通入大量自由氣泡，但這種通氣式微氣泡無法穩定留在近壁區邊界層.微氣泡研究主要集中于減阻效應，氣液兩相流演化非常復雜，尚缺乏完善的減阻機理理論[14-19].Zhu 等[20-21]在前期已完成微氣泡陣列流動減阻性能試驗的研究，提出自穩式電解微氣泡陣列氣膜減阻技術；遴選出較優的微柱孔陣列形狀（直徑d=250 μm，間距比d/l=1.25，徑深比d/h=2），在0.3～1.0 m/s 流速下減阻率可達18%～30%.本文采用大渦模擬(large eddy simulation,LES)，構建平板及駐留微氣泡2 種模型下的近壁區邊界層湍流數值模型.基于湍流數值計算數據，利用本征正交分解(POD)方法進行模態分析，提取各階模態所占的能量和對應的流場結構，開展駐留式微氣泡近壁面湍流流場演化特性、湍流特征參數及減阻機理分析，揭示駐留式微氣泡對湍流邊界層內流場演化特性的影響機制，完善微氣泡減阻機理.

1 數值模擬

1.1 數值求解方法

由于湍流流動為非線性的復雜流動，流體內存在大小不同的渦流和渦旋，大多數工程問題中的流體流動均處于湍流狀態.一般認為即使湍流運動非常復雜，非穩態的連續方程和Navier-Stokes 方程仍然適用于湍流的瞬時運動[22].Navier-Stokes 方程的非守恒形式為

式中：?為通量，ρ為流體密度，t為時間，u、v、w為瞬時速度在x、y 和z 方向的分量，Γ為廣義擴散系數，S 為源相.

目前，湍流數值模擬方法可以分為直接數值模擬方法(DNS)和非直接數值模擬方法，其中非直接數值模擬方法包括雷諾平均法(RANS)及大渦模擬(LES).大渦模擬是介于DNS 與RANS 方法之間的湍流數值模擬方法[23]，基本思想是直接模擬大尺度湍流運動，小尺度湍流表示為亞格子尺度（SGS）湍流模型.為了解析湍流下平板及駐留式微氣泡陣列近壁邊界層區域的精細流場特性，采用大渦模擬作為流場特性分析所使用的湍流模型.

在LES 中，SGS 模型為未分解的湍流運動提供了封閉初始解.根據Smagorinsky[24]的基本SGS 模型，假定SGC 具備以下形式，得到Smagorinsky-Lilly 模型.

式中：τij為亞格子尺度應力；τkk為亞格子尺度應力張量的跡；δij為克羅內克數；μt為亞格子尺度的湍流黏度；Δi為沿i 軸方向的網格尺寸；Cs為Smagorinsky 常數，通常取0.17；為組成的張量.

1.2 平板及駐留微氣泡流動的數值模擬

如圖1 所示為光滑平板及帶有駐留微氣泡陣列表面平板計算域及網格劃分示意圖，計算模型區域大小為L×H=70 mm×10 mm，駐留式微氣泡由微柱孔底部通氣形成.微柱孔參數采用文獻[25]所述，其中微柱孔深度h=125 μm，直徑d=250 μm，間距l=200 μm，計算域內共布置9 個微氣泡陣列.采用多塊結構化網格，對湍流邊界層及微氣泡附近區域網格進行加密（近壁面法向第1 層網格高度y+＜ 1，網格增長率為1.05），整體網格數為4.5×106.

圖1 平板及駐留微氣泡的計算域與網格劃分Fig.1 Computational domain and mesh division of plate and resident microbubbles

為了使邊界層能夠在較短距離內模擬完全發展的湍流結構，采用Sandham 等[26]提出的在邊界層內統一的隨機函數擾動的方法生成湍流邊界層，計算域的邊界條件設置如下.其中計算域入口設置為速度入口，使用用戶自定義函數（UDF）自編譯平均入口速度與內外層擾動速度疊加的流向速度.入口平均速度曲線方程如下所示：

式中：uτ為壁面剪切速度，y 為壁面高度，ν為流體運動黏度.

入、出口流向及展向擾動速度函數如下.

式中：y+為壁面無量綱高度，y+=yuτ/ν ；yp,j為擾動達到峰值時的壁面高度，下標j 為擾動模式序號；ci,j為擾動常數；ωj為擾動頻率，φj為相位移，j=0、1、2、3.

計算域出口設置為壓力出口，上壁面采用自由滑移壁面邊界條件，下壁面采用無滑移壁面邊界條件.

微氣泡水下動態形變過程涉及氣/水兩相流演變，故在數值計算時采用基于體積分數法(volume of fluid,VOF)的氣/水兩相流模型.本文基于大渦模擬(large eddy simulation,LES)模型、VOF 氣/水兩相流模型，構建平板及駐留式微氣泡流動數值模型.在FLUENT 軟件中采用壓力基隱式瞬態求解器，控制方程在離散計算單元中求解，采用SIMPLEC 算法求解壓力-速度耦合，壓力項空間離散采用PRESTO!格式，動量、湍動能及湍流耗散率均使用2 階迎風格式.為了保證計算精度，將非定常計算時間步長設置為0.000 1 s，每個時間步迭代25 次以保證解的收斂，且庫朗數小于0.5.設置計算總時長為0.2 s，使得流動達到穩定完全發展湍流狀態后，提取0.2 s 的數據樣本進行統計分析.

1.3 數值模擬有效性的驗證

湍流邊界層的平均速度分布有3 個典型的區域：線性底層、過渡區和對數區.線性底層與對數區均存在強烈的湍流脈動.在平板邊界層發展為完全湍流后，對整體流場在時間、空間尺度上采取時均統計學處理，提取流向位置x=20δ 時（其中δ 為邊界層厚度，δ=1.45 mm）沿壁面法向上的平均流向速度與相同雷諾數下（Reθ=1 430，θ 為邊界層動量厚度）De Graaff 等[27]的試驗數據進行比較.圖2 給出時均統計后平板湍流邊界層內沿法向的平均速度分布.圖中，＜ ＞表示非定常變量的時間平均結果，U+為無量綱速度，

圖2 平板湍流邊界層內法向平均速度與試驗數據的對比
Fig.2 Comparison of normal average velocity in plate turbulent boundary layer with test data

式中：ρw為水的密度，τ0為壁面切應力.

如圖2 所示，平板湍流邊界層的數值計算結果與試驗值吻合較好，平均誤差小于20%.誤差的產生原因主要有介質不同（空氣與水）、環境因素影響（溫度、氣壓）及仿真模型維度不匹配，導致數值計算結果與試驗數據有一定的差異.

圖3 給出平板湍流邊界層內的各向湍流脈動強度，采用統計學方式提取脈動速度的均方根(RMS)，反映了湍流的強度.圖中，urms+、vrms+分別為流向及法向的無量綱湍流脈動速度的均方根.從圖3 可以看出，平板邊界層的各向湍流脈動強度呈現先升后降的趨勢，在邊界層邊緣逐漸下降至0，峰值均出現在y/ δ=0.022 左右的位置；邊界層內區的脈動強度高于外區，這是由于內區速度梯度大，流速變化較大，使得內區流體受到更強的切應力作用，促使流體出現較強的脈動.將湍流脈動強度與Spalart 等[28-29]的直接模擬(DNS)結果進行比較，平均誤差均小于15%，結果保持一致，驗證了所用數值方法的可靠性.

圖3 LES 與DNS 湍流脈動強度的對比Fig.3 Comparison of turbulent pulsation intensity between LES and DNS

2 POD 方法

POD 方法是將高維數據降維的方法，它的核心是基于數據的協方差矩陣進行特征值分解和特征向量提取的過程.具體來說，對于包含 n 個樣本或快照的數據集合，可以將其表示為 m×n的矩陣A，其中 m 為每個樣本的維度.計算 A的協方差矩陣 C=ATA，對其進行特征值分解.通過特征值和特征向量的排序和截斷，可以將數據集合的信息壓縮到低維的子空間中.這個子空間中的基向量由協方差矩陣的前r 個特征向量構成，其中 r 為要保留的主要特征向量的數量.采用這些基向量來表示原始數據集合中的每個樣本，得到其在低維子空間中的表示.由于本征向量是正交的，在低維子空間中表示的數據具有良好的可解釋性和可視化性.

通過所提取的各瞬時時刻流場信息（速度）集計算出時均流場信息，將瞬時流場信息減去時均流場信息，獲得每個時刻的脈動速度矩陣.采用POD 方法，將脈動速度矩陣分解為正交模態和時間系數，如下所示：

式中：u′(X,t)、u(X,t)和(X,t)分別為某時刻下流場脈動速度矢量場矩陣、瞬時速度矢量場矩陣和時均速度矢量場矩陣；ak(t)為該流場時間系數；φk(X)為分解得到的第k 階模態；N 為提取的流場時刻樣本數，N=200，保證了POD 方法湍流特征統計量的穩定性.

對該流場脈動速度矢量場矩陣分解得到的對應正交基的能量排序，可得第k 階POD 模態的能量與總流場能量的占比，如下所示：

式中：λk為第k 階模態對應的特征值.

3 計算結果及POD 分析

3.1 平板流場及駐留微氣泡流場的演化特性

雖然傳統通氣式微氣泡減阻方式能夠達到約30%的減阻率，但存在減阻效果不穩定及需要不斷通氣的問題；駐留微氣泡陣列兼具良好的穩定性與減阻性能.圖4 給出流場演化過程中平板及駐留微氣泡陣列表面在x=15δ～25δ 位置的壁面切應力τ 分布曲線，其中無量綱時刻t+=t/t1（t1為流場演化的時間0.2 s，t 為流場演化過程中的各個時刻）.如圖4 所示，駐留微氣泡陣列表面湍流邊界層平均壁面切應力約為1.23×10-3N，相較于平板減少了約13.7%，驗證了駐留微氣泡陣列具備一定的減阻性能.在微氣泡氣/水界面動態形變及湍流特性的影響下，表面壁面切應力呈現動態振蕩特性，相較于平板在t+≈ 0.52 時刻下的突變，微氣泡陣列表面的切應力表現得更穩定.

圖4 平板及微氣泡陣列表面壁面切應力分布Fig.4 Wall shear stress distribution of plates and surface of microbubble arrays

以往對近壁面湍流相干結構的研究表明，在一對正、負的流向渦旋之間存在低速條帶[30]，如圖5 所示.流向渦旋上沖側帶離壁面的低速流體稱為上拋，流向渦旋下沖側帶向壁面的高速流體稱為下掃，這些向上和向下的流體運動被稱為湍流相干結構猝發.流向渦結構內高速流體的下掃產生的近壁面雷諾切應力是近壁面湍流摩擦阻力的主要來源.

圖5 上拋、下掃猝發[30]Fig.5 Upcast,down-sweep burst[30]

圖6(a)～(d)給出不同時刻平板湍流邊界層流向速度場u 的演化云圖.如圖5 所示，完全發展的湍流平板邊界層內速度分布不再呈現時均流場下穩定的邊界層，而是形成多個低速流體突起與高速流體下凹，平板湍流邊界層存在低速回流渦結構.如圖6(e)所示為該低速區速度矢量圖.可以看出，該區域呈現低速流體上拋與高速流體下掃的“猝發”現象，形成湍流相干結構.從圖6 還可以看出，隨著時間的遷移，速度分布突起從上游逐漸遷移至下游，當t+≈ 0.5 時，流向渦結構遷移至x=15δ～25δ 位置，與平板切應力突變節點相對應.

圖6 不同時刻平板湍流邊界層流向速度場Fig.6 Flow velocity contour of plate turbulent boundary layer at different time

圖7(a)～(d)給出駐留式微氣泡陣列表面不同時刻流向速度u 的演化云圖.可以看出，駐留式微氣泡陣列的存在使得整體速度場更均勻，邊界層內流向渦在流經微氣泡后消失，表明微氣泡陣列可以破壞湍流相干結構.如圖7(e)所示為微氣泡陣列表面速度矢量.可以看出，駐留微氣泡動態形變下其柔性氣/水界面處存在豐富的渦系結構，微氣泡前部與后部存在對稱、旋向相反的渦旋，前部形成的逆時針旋轉渦旋表現為高速流體上拋、低速流體下掃，抑制“猝發”.微氣泡后部形成的渦旋方向雖然與相干結構一致，但其直徑較平板流向渦減小了約80%.從圖7(e)還可以看出，氣/水界面動態形變使其曲率不斷變化，微氣泡表面微曲率導致的康達效應減弱，微氣泡后部形成流動分離，在下一個微氣泡前部形成再附著.如此產生的渦系結構極大程度地破壞了湍流邊界層分布的連續性，減少了微氣泡陣列邊界層內高摩擦阻力的產生.

圖7 不同時刻微氣泡陣列表面流向速度云圖Fig.7 Flow velocity contour of surface of microbubble array at different time

圖8 給出不同時刻平板及駐留微氣泡陣列表面在湍流邊界層內(x=20δ，y+=10)流向及法向的無量綱速度u+、v+演化曲線，據此分析微氣泡陣列對湍流相干結構“猝發”的量化影響機制.如圖8所示，平板湍流邊界層內流向速度及法向速度的幅值均大于微氣泡陣列表面，存在明顯的周期性演化.該速度演變周期為湍流相干結構“猝發”周期.通過對法向速度演化實施傅里葉變化(FFT)，得到頻譜圖，如圖9 所示.圖中，vma為平板及微氣泡陣列表面法向速度的幅度.從圖9 可以看出，平板法向速度主頻率(8.77 Hz)與主頻所在幅值(59.87 m/s)均大于駐留微氣泡陣列表面主頻率(3.17 Hz)與幅值(16.78 m/s)，說明微氣泡陣列的存在將“猝發”頻率有效降低了63.85%，達到了約13.7%的減阻效果.

圖8 不同時刻平板及微氣泡陣列表面速度演化Fig.8 Velocity evolution of plates and surface of microbubble arrays at different time

圖9 平板及微氣泡陣列表面法向速度頻譜圖Fig.9 Spectrum of normal velocity on plate and surface of microbubble array

如圖10 所示為平板與微氣泡陣列表面邊界層近壁面相同區域內統計得到的平均湍動能TKE 演化曲線.可以看出，微氣泡陣列結構使得近壁面湍動能較平板減少(微氣泡陣列氣膜平板表面平均湍動能為0.093 m2/s2，平板表面平均湍動能為0.165 m2/s2)，減少了湍流與近壁面之間的能量交換，降低了流體與壁面之間的相互作用.

圖10 平板及微氣泡陣列表面湍動能演化Fig.10 Turbulence kinetic energy evolution of plate and surface of microbubble array

3.2 流場演化POD 分析

采用POD 方法，對應用大渦模擬方法得到的湍流邊界層流場特征結構進行降階處理.選取流場范圍為長度L1=40 mm，高度H1=6 mm，流場演化時間t1=0.2 s，各個瞬時時刻間隔Δt1=0.001 s，將時間劃分為200 個樣本.流向速度u 與法向速度v 的樣本空間大小為1 001×601，生成201 個快照向量矩陣，從而得到201 個POD 基向量.通過對流場二維速度矢量場進行POD 分解，求解得到流場速度的旋度（即渦量），分析湍流邊界層內擬序結構的演變規律.

圖11 給出基于POD 方法分解得到的光滑平板與微氣泡陣列表面湍流邊界層內第0 階模態（平均流場）云圖.圖中，ωa為渦量.從圖11 可以看出，采用POD 分解得到的平均流場云圖與時均化流場邊界層結構高度一致，流場穩定性較強.在光滑平板邊界層內，總體渦量強度相較于微氣泡表面高.在微氣泡陣列區域，渦量強度較高，表明微氣泡的存在破壞了湍流邊界層內原有的結構.

圖11 0 階模態/平均流場Fig.11 Zero order/average flow field

圖12 給出二維平板速度矢量場經POD 分解后得到的前20 階模態總能量累積占比Pc分布曲線.可以看出，平板邊界層達到總體湍流能量的90%需要11 階的模態，而微氣泡陣列表面只需要前7 階模態，這說明微氣泡的存在使得平板邊界層內的湍流流場結構更穩定，總體減少了湍流擬序結構.圖13 給出單獨各階模態的能量占比Pk情況，其變化規律隨階數的增加而減少，最終趨于0.將平板二維速度矢量場分解至10 階模態，可以完整重構原有流場.

圖12 各階總能量累積分布曲線Fig.12 Cumulative distribution curve of total energy of each order

圖13 各階能量占比變化曲線Fig.13 Variation curve of energy proportion for each order

圖14、15 分別給出完全發展的平板湍流邊界層與微氣泡陣列氣膜邊界層二維速度矢量場前10 階模態分解云圖.可以看出，隨著模態階數的增加，平板湍流邊界層內渦旋結構逐漸失穩，渦量強度增大，到第7 階模態出現典型的湍流相干擬序結構，渦團結構呈現正負交替對稱分布，遍布整個平板邊界層，與3.1 節所述流向渦的“猝發”過程一致，形成高摩擦阻力.對比微氣泡陣列氣膜下各階模態的分解云圖可知，微氣泡的存在對整體平板擬序結構具備強烈的破壞性，微氣泡陣列氣膜表面渦旋結構形態更扁平，湍流特征結構在壁面法向陣列表面的擬序結構較光滑平板明顯減少，相應受到抑制.對比圖14(e)、15(e)的第9 階模態可見，微氣泡近壁區的切應力降低，表現出良好的減摩特性.

圖14 平板湍流邊界層的模態云圖Fig.14 Modal nephograms of turbulent boundary layer on plate

圖15 微氣泡陣列表面湍流邊界層的模態云圖Fig.15 Modal nephograms of turbulent boundary layer on surface of microbubble array

圖16 給出典型階數下微氣泡陣列氣膜的模態局部云圖.可以看出，微氣泡陣列局部大尺度渦旋結構顯著減少，加強了湍流區域內的小尺度結構，使得流場內的湍流動能分布更均勻，抑制了湍流擬序結構的發展，在邊界層內形成良好的減阻效果.

圖16 局部微氣泡陣列表面湍流邊界層的模態云圖Fig.16 Modal nephograms of turbulent boundary layers on surface of local microbubble arrays

4 結 論

(1)微氣泡陣列（直徑d=250 μm、間距比d/l=1.25 及徑深比d/h=2）的存在使得近壁區切應力減小約13.7%，流向渦尺度降低約80%，削弱了動量交換，降低了湍動能與動量耗散.

(2)流向渦的發展、遷移帶動邊界層內低速流體上拋、高速流體下掃，形成“猝發”現象，是平板湍流邊界層內高摩擦阻力的本質原因.

(3)駐留微氣泡的存在可以破壞湍流相干結構，使得“猝發”頻率由8.77 Hz 降低至3.17 Hz，脈動速度幅值隨之降低.

(4)利用POD 方法，能夠有效地提取湍流擬序結構，微氣泡的存在使得平板表面渦旋結構形態更扁平，加強了湍流區域內的小尺度結構，使得流場內的湍流動能分布更均勻，抑制了相干結構的發展.將駐留微氣泡應用于水下航行體以提升表面減阻性能是下一階段研究的重點.