分層流體中冰脊拖曳力的數值模擬研究

祖永恒，盧 鵬，王慶凱，李志軍，吳 巖，李 博

（大連理工大學 港口海岸及近海工程國家重點實驗室，遼寧 大連 116024）

1 研究背景

流體分層是自然界普遍存在的一種現象，包括大氣分層和水體分層。流體分層主要受其溫度和密度的控制，在以密度為主導的分層環境中，物體在流體中的相對運動會呈現出不同于均勻流的動力特性。例如，河口附近由于海水和淡水的密度差異會形成一個切變鋒面，鋒面上層是淡水層，下層是海水層，這對通行的船舶以及寒區河冰的運動都有重要影響［1-3］。冰蓋或冰塞等對上層流體的擾動會在流體分界面激起內波，內波同時又會改變冰蓋或冰塞受到的拖曳力，從而進一步影響冰塞的堆積演變。王軍等研究了均勻流體條件下直槽中冰塞的形成機理，并給出了試驗條件下不同冰塞形成機制的臨界弗勞德數Frc范圍在0.12～0.13［4］。寒區水庫解凍初期的水溫分布差異也會引起表層的浮力流動現象，對水庫的冰情、下泄水溫造成影響［5］。李嘉等使用二維水庫水溫模型分析了水庫的分層流場和溫度場的耦合規律［6］；鄭鐵鋼等采用量綱分析方法研究了水溫分層對水庫下泄水溫的影響［7］。但是寒區水庫中冰凌冰塞的產生也會改變分層流體的水動力情況，進而也會對水庫下游取水的水溫、水質造成影響，針對這方面的相關研究有待深入。

海洋中的流體分層現象也較為普遍，特別是近年來北極夏季海冰加速融化，導致海洋鹽躍層變淺甚至在海冰邊緣區和冰脊深度達到同一量級，此時海冰運動激發界面內波將會對冰脊拖曳力造成不可忽略的影響［8-10］，甚至發生死水現象［11］。根據海冰拖曳系數參數化的思想，冰-水總拖曳力包括摩拖曳力和形拖曳力兩部分。其中，摩拖曳力描述的是由海冰表面均勻分布的小凸起物引起的剪切力，形拖曳力描述的是由海冰表面非均勻分布的較大凸起物（如冰側、冰脊）引起的流場變化而導致的水平方向壓力差［12］，不同類型的冰面粗糙物由截斷高度來劃分［13］。冰脊形拖曳力對冰-水總拖曳系數的貢獻可由下式定義：

式中：A 為冰密集度；hk/Lr為冰脊密度；即冰脊的入水深度hk和冰脊間距Lr之比；Cr為單個冰脊的局地拖曳系數。若A 和hk/Lr一定，Cr直接決定了冰脊對冰-水界面拖曳系數的貢獻。在冰脊分布密度較大時，冰脊形拖曳力對冰-水總拖曳系數的貢獻不可忽略［14］。因此，研究單個冰脊局地拖曳系數的變化規律對于冰-水界面的動量交換研究具有重要意義。

真實的流體分層存在連續分層和多層等復雜情況，可以采用雙層流體對該問題進行簡化。雙層流體的物理實驗具有較強可操作性，可以精確再現流體分層情況，并反映出不同流態下分層流體的變化特征。Pite 利用該方法對冰脊拖曳力問題進行了研究，并結合Baines 對內波流場的描述，根據內弗勞德數（無量綱拖曳速度）和冰脊入水深度的變化，將雙層流流況分為亞臨界區、跨臨界區（下風波區）和超臨界區［15-16］。Jammel 利用基于歐拉方程的有限差分算法對Pite 的分層流體實驗進行了數值模擬［17］；Mortikov 利用浸沒邊界法對雙層流體的邊界進行了更精細化的處理［18］。這些研究多是對雙層流產生的內波和內波拖曳力進行定性分析，并未建立雙層流中冰脊形態和冰-水拖曳系數的參數化關系；而且已有數值模擬研究中很少考慮湍流模型的應用，對冰脊后的流場壓強也缺少分析。隨著計算流體力學的發展，通過并行計算已經大大提高了湍流模型的計算效率。本文采用基于有限體積法的RNG k-ε湍流模型和VOF 界面追蹤方法對該問題進行數值模擬研究。通過對數值模擬流場的處理，與物理實驗的結果進行了對比，分析了動壓強對內波增阻的影響，揭示了冰脊局地拖曳系數隨著冰脊入水深度和流場弗勞德數變化的一般規律。

2 分層流試驗研究

數值模擬計算域的設置與先期進行的實驗室物理實驗在垂向尺寸上完全一致。物理實驗是在大連理工大學海岸及近海工程國家重點實驗室的粒子圖像測速水槽完成（見圖1），該實驗的詳細論述見文獻［19］。水槽長450 cm，寬23 cm，深45 cm，上層淡水深15 cm，下層鹽水深20 cm。冰脊模型的材料為有機玻璃，冰脊底角設計為45°。拖車由電機牽引在水槽上方的軌道上運行，可帶動冰脊以1、3、6、7、8、9、10、12、15、18、24、30 cm/s 等不同的速度在水槽內勻速運動。冰脊入水深度可調，實驗設置了4、6、8、10 cm 4 種入水深度。冰脊模型在水平方向上連接拉壓傳感器，再固定在拖車上，能夠測量模型在水平方向受到的流體拖曳力。

圖1 物理實驗水槽

水槽中產生的內波最大傳播速度不會超過其相速度。本研究中的內波相速度約為14 cm/s，因此數值模擬中計算域長度分別向上下游延長實際水槽長度的一倍，數值計算時間不超過40 s 以消除側邊界反射的影響。對于流體的物理參數，數值模擬和物理實驗的取值保持一致，淡水的密度為998.2 kg/m3，動力黏性系數為1.0003×10-3kg/（m·s）；海水的密度為1025 kg/m3，動力黏性系數為1.12×10-3kg/（m·s）。考慮到在水槽寬度方向上流速變化較小，將數值模擬簡化為二維問題，冰脊以恒定的速度和入水深度在水面處沿水平方向運動，流體初始速度為0，和物理實驗的組次一致，見圖2。

圖2 數值模擬水槽示意圖（單位：cm）

本研究模擬雙層流中內波的產生，遵循重力相似準則。空氣、上層流體密度差與上下層流體的密度差的比值為（ρ1-ρa）/（ρ2-ρ1）≈37，產生的表面波的波高和內波的波高之比也為1/37，因此表面波的影響忽略不計［20］。弗勞德數Fr 是重力相似準則中判定流況動力相似的無量綱參數，定義為冰脊拖曳速度和線性內波相速度的比值：

式中：V 為冰脊拖曳速度；C 為線性內波相速度。模型比尺應滿足：

內波相速度采用線性非色散的長波相速度［18］：

式中ρ0、h0為模型中的特征密度和特征深度。Pite 統計的現場觀測數據顯示，當冰脊運動速度在0.2 m/s、入水深度為10 m 時，弗勞德數的變化范圍是0.12 ～0.69，冰脊入水深度和上層流體深度之比的變化范圍是0.12 ～1.38；但在后來的觀測資料中，弗勞德數的變化范圍可以達到0.2 ～3.8［16，22］。數值模擬和物理實驗中弗勞德數的變化范圍是0.07 ～2.14，冰脊入水深度和上層流體深度之比的變化范圍是0.27 ～0.67，保證了本研究結果能夠覆蓋現場不同的內波形態，包括亞臨界區、跨臨界區和超臨界區。

在數值模擬的計算域中（見圖2），上邊界DK 為壓強入口，入口壓強為0（參考壓力為大氣壓）；外邊界DG、GH、HK 都為壁面，速度為0；其中冰脊ABC 也是壁面，同時也是運動的內邊界，運動速度為冰脊的拖曳速度。

計算域用三角形網格和四邊形網格進行構建。圖3是網格劃分示意圖，其中冰脊附近的區域為三角形非結構網格，其他區域是四邊形結構網格。包含冰脊部分的網格設定為動網格，冰脊邊界以恒定速度向左運動，下方為靜止網格。動網格區域和靜網格區域通過一對interface 界面實現數據交換。

表1 網格無關性驗證結果

圖3 網格劃分示意圖（紅色為水-氣界面，藍色為密度分界面）

數值模擬中的冰脊入水深度最大為10 cm，非結構網格也最復雜，故網格大小的無關性驗證選取T=10 cm、V=15 cm/s 組次，檢驗結果見表1。當網格步長為0.8 cm ～1.2 cm 時，計算結果已趨于穩定，相對誤差不超過3%。故選取空間步長1 cm。

流場計算的控制方程采用基于雷諾應力平均方程的RNG k-ε湍流模型，k-ε兩方程如下：

在Fluent 軟件中，Cμ=0.0845，C1ε=1.42，C2ε=1.68，αε=αk=1.39。控制方程在時間離散方式上采用一階隱式離散，空間離散方式上采用二階迎風離散。求解算法是在SIMPLE 算法（壓力耦合方程組半隱式方法）基礎上修正的PISO 算法（隱式算子分割算法）。冰脊的運動采用動網格模塊，使用鋪層（layering）方案進行求解。對于界面的追蹤，選用基于Youngs 算法Geo-Reconstruct 方案計算網格邊界的流體體積量。

3 計算結果和分析

為了方便分析，引入無量綱參數B，定義為冰脊入水深度T和上層流體深度h1的比值：

對于單個冰脊的局地拖曳系數Cr，采用阻力系數的定義形式：

式中：F 為冰脊受到的流場拖曳力；ρ1為上層流體密度；A 為冰脊在垂直于運行方向的投影面積；V為冰脊的運動速度。

3.1 冰脊拖曳力隨速度的變化冰脊拖曳力F 由冰脊表面應力（包括黏性應力和壓強應力）在水平方向上的積分得到，F 隨冰脊運動速度V 的變化情況以及與物理實驗結果的對比如圖4所示。

圖4 數值模擬冰脊拖曳力和物理實驗對比結果

圖5 數值模擬拖曳系數隨弗勞德數的變化和物理實驗的對比結果

由圖4可以看出來，數值模擬和物理實驗在Fr＜1 時均表現出了非線性，拖曳力隨速度的增加先增大后減小。當Fr＞1 時，拖曳力隨著速度的增大而增大。這一現象在拖曳系數的變化規律中更加明顯，而且當Fr＞1 后，拖曳系數趨近于1，和單層流拖曳系數的變化規律相似，見圖5中數值模擬和物理實驗的對比結果。

在圖5（b）中，物理實驗得到的拖曳系數的非線性變化并不明顯，這是由于物理實驗條件的限制造成的。當速度為1 cm/s 時，物理實驗中的拖曳力值甚至不到10-3N，由力值計算得到的拖曳力系數也會產生誤差，導致了物理實驗在小速度區間的拖曳系數變化規律不太明顯，但此時數值模擬中拖曳系數的非線性變化則較為顯著。當Fr＜1 時，拖曳系數隨速度先增大后減小，峰值隨著入水深度的變化而變化。當Fr＞1 后，拖曳系數不再隨速度的變化而變化，而是趨近于單層流體的冰脊拖曳系數。

需要注意的是，圖4中的數值模擬的冰脊拖曳力結果相比物理實驗存在一定差異，主要是由兩種方法中流體密度剖面的差異造成的，密度剖面的差異導致了浮頻率N 的變化，見圖6數值模擬中流體分層條件和物理實驗的對比。在流體分界面，數值模擬可以直接定義一個密度斷面，斷面上浮頻率趨近于無限大。但物理實驗需要在一個深度區間內完成密度的連續變化，即在約7 cm 的區間內由淡水過渡到鹽水，圖6（b）物理實驗中的浮頻率最高達到3 rad/s。

對于密度線性變化的連續分層流體，當h=35 cm 時，其內波相速度C：

而本研究采用兩層流體線性假設下的內波相速度C：

圖6 數值模擬中流體分層條件和物理實驗的對比

由上面兩式可知相速度C 和密度梯度相關，數值模擬中分界面的密度梯度大于物理實驗值；而物理實驗中在分界面附近屬于連續分層，在整體上又屬于雙層流情況，因此物理實驗的內波最大相速度應小于數模，介于8.1 cm/s 和14.4 cm/s 之間，并偏于14.4 cm/s。所以圖4中物理實驗拖曳力的非線性區間比數值模擬值提前，該區間拖曳力的峰值也比后者提前，造成了數值模擬和物理實驗出現差異。

當Fr＞1 時，雙層流中冰脊拖曳力受內波影響較小，內波拖曳力幾乎為零，拖曳力的變化規律和單層相似，此時數值模擬和物理實驗吻合良好。單層流冰脊拖曳力數據可以參考吳巖等人的研究成果［23］。另外，圖7（a）給出了冰脊入水深度T=8 cm 時，在雙層流中的拖曳力與在單層流中拖曳力數值模擬結果的對比，從圖中可以看出相同水深情況下，雙層流的阻力情況相對單層流復雜一些。以Fr=1為分界點，前面區間的阻力差異較大，這是由于內波的生成，對冰脊后的壓力場產生了較大的影響，本文后面章節還會詳細分析。當Fr＞1 后，內波對冰脊阻力的影響幾乎可以忽略了，單、雙層流中冰脊拖曳系數都在1 左右，見圖7（b）。

圖7 T=8 cm 時，數值模擬的冰脊拖曳力和拖曳系數在單層和雙層流中的變化結果

3.2 冰脊拖曳力隨入水深度的變化雙層流中冰脊拖曳力、拖曳系數隨無量綱入水深度B 的變化規律見圖8所示。

當Fr 一定時，拖曳力隨著入水深度的增加而增大，接近于線性變化。從拖曳系數來看，當Fr＞0.57 時，拖曳系數的變化較小，接近于1，如圖8中Fr=0.86、Fr=1.07、Fr=2.14 時的情形。當Fr≤0.57 時，拖曳系數隨著入水深度的增加而顯著增大，而且拖曳速度越小，拖曳系數變化幅度越大，如圖8中Fr=0.42、Fr=0.57 的情形。

由拖曳系數的定義可知，拖曳系數和冰脊運動速度的平方呈反比。因此當速度較小的時候，拖曳系數對拖曳力的變化更加敏感。拖曳力的非線性變化出現在弗勞德數較小的區間，此時拖曳系數隨Fr 的變化趨勢更為顯著（圖7）。因此圖8中Fr 較小時拖曳系數的變化幅度要大于Fr 較大時的變化幅度。

3.3 內波形態的變化規律圖9給出了入水深度T=6 cm，速度V=15 cm/s 時冰脊運動對應的典型內波形態。在拖曳系數隨Fr 變化的非線性區間，內波波速大于冰脊的運動速度，因此內波的擾動可以向上下游同時傳播。此時內波向下游傳播到一定的距離便完全耗散，而內波波谷向上游延伸，并且波面變化幅度較小，近似平緩的曲線。在下游形成的內波比較穩定，波峰相對冰脊的位置都保持不變。

圖8 冰脊拖曳系數隨無量綱入水深度的變化結果

圖9 冰脊入水深度T=6 cm，速度V=15 cm/s 數值模擬的內波

圖10 數值模擬內波無量綱化波高和物理實驗對比結果

圖11 冰脊入水深度T=8 cm，速度V=7 cm/s 數值模擬內波形態和物理實驗的對比結果（黑色虛線為初始分界面）

隨著弗勞德數的增長，冰脊上游的擾動消失，僅冰脊下方出現波谷狀態。此時，下游內波仍處于比較穩定的狀態。而當Fr＞1 時，由于冰脊的拖曳速度超過了內波的相速度，內波擾動僅向冰脊下游傳播，而且第一個波峰逐漸變得平坦，并逐漸遠離冰脊。當Fr=2.1（V=30 cm/s）時，內波形態變得更加雜亂，波峰不斷遠離冰脊，最終擾動逐漸消失，流態過渡到超臨界區域。

其他入水深度條件下產生的內波形態和T=6 cm 的情形一樣，但波峰高度有所不同。為進行數值模擬和物理實驗的比較，把內波最低點與最高點的高度差定義為內波波高H，H/h1為無量綱化波高，圖10是不同入水深度情況下無量綱波高的對比。

圖10中數值模擬的波高比物理實驗結果偏大，但趨勢一致。在跨臨界狀態，波高隨著速度的增大而增大；在超過臨界狀態后，波高開始趨于平穩。數值模擬與物理實驗的差異還是由于上文中提到的密度分層的不同，它導致了圖11中物理實驗的摻混現象比較明顯，消耗了更多的能量。而數值模擬中的波形則比較穩定，因此內波波高相對較大。

3.4 內波對冰脊拖曳力的影響由圖7中的對比可知，在拖曳力變化的非線性區間，內波對冰脊拖曳力產生了重要影響，下面從垂直于冰脊運行方向和平行于冰脊運行方向兩個方面對該問題進行分析。首先對于垂直方向，統計了非線性區間的流體界面垂直偏移的最大值D，該距離隨弗勞德數的變化規律見圖12。從圖中可知，偏移距離D 和拖曳力F 在跨臨界區域相關性較高，相關系數R2＞0.94。分界面位移的變化可以反映出流體勢能的變化，同時，由于上層流體通道束窄，流速也會增加。同一流線在該位置勢能增加，動能增加，而且下層流體對上層流體的剪切力做負功，流線的壓強勢能則會降低，在分界面位移最高處形成一個相對低壓區。

圖12 數值模擬分界面位移和冰脊拖曳力關聯圖（藍線是拖曳力，黑線是位移）

對于平行冰脊運動方向，冰脊拖曳力主要是由冰脊前后壓強場的變化引起的。因此對冰脊附近的壓強場進行后處理，將流場總壓強減去靜壓，得到流場的動壓強。這樣更容易反映出冰脊前后壓強的變化，結果見圖13。

圖13是不同運動速度冰脊附近的壓強場和流場。可以看到雙層流體中不僅會在冰脊后形成一個尾流漩渦場（和單層相同），而且還會在內波擾動的波峰、波谷處形成漩渦場。每一個漩渦場都是一個相對低壓區，渦動中心的壓強最低，這也是冰脊產生拖曳力的原因。在單層流中，冰脊拖曳力只是受尾流漩渦場一個低壓區的影響，而且隨著速度的增大，漩渦中心的相對壓強降低。然而在雙層流中，冰脊拖曳力要受尾流漩渦場和冰脊后第一個內波波峰漩渦場的雙重影響。尤其當拖曳速度小于內波相速度時，內波擾動向冰脊上下游同時傳播，此時內波的形態相對冰脊比較穩定，跟隨冰脊的運動而運動。波峰的漩渦場和尾流漩渦場共同作用，對冰脊拖曳力產生影響。

圖13 數值模擬分層流體冰脊附近動壓強分布云圖（黑色實線為內波波面）

隨著速度的增大，圖13中內波波峰不斷遠離冰脊，對冰脊尾流場的影響也越來越小，冰脊拖曳力也逐漸趨近于單層流的情況。在Fr=0.71（V=10 cm/s）時，內波波峰的漩渦區對冰脊尾流的漩渦區影響最大，內波波峰的高度也達到最高（圖12（b））。此時尾流漩渦區不能充分發展，受到內波波峰漩渦的“擠壓”效果，漩渦的相對壓強也最低，因此出現了冰脊拖曳力的極值。其他入水深度條件的冰脊拖曳力變化規律和T=6 cm 的情況類似。

4 冰脊底角的影響

除了流速和入水深度會對冰脊的拖曳力產生影響外，冰脊形狀也會造成一定的差異。利用數值模擬可以探究冰脊拖曳力隨著不同冰脊形狀或者分層水深的變化規律，但由于篇幅限制，只給出冰脊形狀變化的定性結果。

由圖14可知，當冰脊底角在30°到60°范圍內，拖曳力的變化規律基本相似。相同速度和入水深度下，拖曳力角度的增大而變大。冰脊拖曳系數隨著角度的增加而增加，在非線性區間向線性區間的轉換過程中，拖曳系數隨角度的變化稍微平緩，見圖15中Fr=0.71 的情況，其他弗勞德數區間的變化趨勢基本一致。內波形態的演變對冰脊底角的變化并不敏感，主要還是受弗勞德數的影響。

圖14 T=6 cm，不同冰脊底角的拖曳力隨弗勞德數的變化

圖15 T=6 cm，拖曳系數隨角度的變化：Fr≤0.71 為非線性區間、Fr＞0.71 為單調區間

5 結論

雙層流環境下冰脊拖曳力隨著冰脊運動速度和冰脊入水深度的變化而變化，對應的冰脊拖曳系數受弗勞德數和冰脊的無量綱入水深度的影響。數值模擬得到的冰脊拖曳力和拖曳系數與物理實驗結果基本吻合，在跨臨界區域內都表現出冰脊拖曳力和拖曳系數的非線性變化規律。數值模擬產生的內波形狀和物理實驗結果一致，但由于分界面密度分布的差異，物理實驗中存在更多的摻混現象，消耗了更多的內波能量。

當入水深度保持不變時，冰脊拖曳力隨著拖曳速度的增加有先增加后減小、然后再增加的趨勢。冰脊拖曳系數隨弗勞德數變化的規律更加明顯，非線性變化峰值出現的位置也較拖曳力提前。當弗勞德數Fr＞1 時，雙層流的冰脊局地拖曳系數穩定于1，與單層流的冰脊局地拖曳系數差異不大。對非線性區間分界面位移和冰脊附近壓強場的分析發現，當位移達到最大時，拖曳力達到一個峰值。而且內波波峰和流場壓強相互耦合，隨著分界面垂直偏移距離的增加，內波波峰附近的負壓區域增大，在波峰位置不隨時間變化的情形下，會影響冰脊后的尾流負壓區，使拖曳力出現峰值。當冰脊運動速度保持不變時，在Fr≤0.57 的較小范圍內，冰脊拖曳系數隨入水深度的增加而增加；在Fr＞0.57 的較大范圍內，冰脊拖曳系數基本不隨入水深度的增加而變化，最終趨近于1。

單個冰脊的局地拖曳系數對冰-水總拖曳系數的確定具有重要作用，結合流體分層的物理實驗研究，本文利用數值模擬方法探討了雙層流中冰脊拖曳系數隨冰脊入水深度和弗勞德數變化關系，為分層流體冰-水拖曳系數的參數化奠定了基礎。除了冰脊入水深度和運動速度，冰脊形態和分層水深的變化也會對冰脊局地拖曳系數產生影響，需要進一步的試驗研究和理論分析，從而發展完整的冰脊拖曳系數參數化方案。