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

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

（大連理工大學(xué) 港口海岸及近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024）

1 研究背景

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

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

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

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

2 分層流試驗(yàn)研究

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

圖1 物理實(shí)驗(yàn)水槽

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

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

本研究模擬雙層流中內(nèi)波的產(chǎn)生，遵循重力相似準(zhǔn)則。空氣、上層流體密度差與上下層流體的密度差的比值為（ρ1-ρa(bǔ)）/（ρ2-ρ1）≈37，產(chǎn)生的表面波的波高和內(nèi)波的波高之比也為1/37，因此表面波的影響忽略不計(jì)［20］。弗勞德數(shù)Fr 是重力相似準(zhǔn)則中判定流況動(dòng)力相似的無(wú)量綱參數(shù)，定義為冰脊拖曳速度和線(xiàn)性?xún)?nèi)波相速度的比值：

式中：V 為冰脊拖曳速度；C 為線(xiàn)性?xún)?nèi)波相速度。模型比尺應(yīng)滿(mǎn)足：

內(nèi)波相速度采用線(xiàn)性非色散的長(zhǎng)波相速度［18］：

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

在數(shù)值模擬的計(jì)算域中（見(jiàn)圖2），上邊界DK 為壓強(qiáng)入口，入口壓強(qiáng)為0（參考?jí)毫榇髿鈮海煌膺吔鏒G、GH、HK 都為壁面，速度為0；其中冰脊ABC 也是壁面，同時(shí)也是運(yùn)動(dòng)的內(nèi)邊界，運(yùn)動(dòng)速度為冰脊的拖曳速度。

計(jì)算域用三角形網(wǎng)格和四邊形網(wǎng)格進(jìn)行構(gòu)建。圖3是網(wǎng)格劃分示意圖，其中冰脊附近的區(qū)域?yàn)槿切畏墙Y(jié)構(gòu)網(wǎng)格，其他區(qū)域是四邊形結(jié)構(gòu)網(wǎng)格。包含冰脊部分的網(wǎng)格設(shè)定為動(dòng)網(wǎng)格，冰脊邊界以恒定速度向左運(yùn)動(dòng)，下方為靜止網(wǎng)格。動(dòng)網(wǎng)格區(qū)域和靜網(wǎng)格區(qū)域通過(guò)一對(duì)interface 界面實(shí)現(xiàn)數(shù)據(jù)交換。

表1 網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果

圖3 網(wǎng)格劃分示意圖（紅色為水-氣界面，藍(lán)色為密度分界面）

數(shù)值模擬中的冰脊入水深度最大為10 cm，非結(jié)構(gòu)網(wǎng)格也最復(fù)雜，故網(wǎng)格大小的無(wú)關(guān)性驗(yàn)證選取T=10 cm、V=15 cm/s 組次，檢驗(yàn)結(jié)果見(jiàn)表1。當(dāng)網(wǎng)格步長(zhǎng)為0.8 cm ～1.2 cm 時(shí)，計(jì)算結(jié)果已趨于穩(wěn)定，相對(duì)誤差不超過(guò)3%。故選取空間步長(zhǎng)1 cm。

流場(chǎng)計(jì)算的控制方程采用基于雷諾應(yīng)力平均方程的RNG k-ε湍流模型，k-ε兩方程如下：

在Fluent 軟件中，Cμ=0.0845，C1ε=1.42，C2ε=1.68，αε=αk=1.39?？刂品匠淘跁r(shí)間離散方式上采用一階隱式離散，空間離散方式上采用二階迎風(fēng)離散。求解算法是在SIMPLE 算法（壓力耦合方程組半隱式方法）基礎(chǔ)上修正的PISO 算法（隱式算子分割算法）。冰脊的運(yùn)動(dòng)采用動(dòng)網(wǎng)格模塊，使用鋪層（layering）方案進(jìn)行求解。對(duì)于界面的追蹤，選用基于Youngs 算法Geo-Reconstruct 方案計(jì)算網(wǎng)格邊界的流體體積量。

3 計(jì)算結(jié)果和分析

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

對(duì)于單個(gè)冰脊的局地拖曳系數(shù)Cr，采用阻力系數(shù)的定義形式：

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

3.1 冰脊拖曳力隨速度的變化冰脊拖曳力F 由冰脊表面應(yīng)力（包括黏性應(yīng)力和壓強(qiáng)應(yīng)力）在水平方向上的積分得到，F(xiàn) 隨冰脊運(yùn)動(dòng)速度V 的變化情況以及與物理實(shí)驗(yàn)結(jié)果的對(duì)比如圖4所示。

圖4 數(shù)值模擬冰脊拖曳力和物理實(shí)驗(yàn)對(duì)比結(jié)果

圖5 數(shù)值模擬拖曳系數(shù)隨弗勞德數(shù)的變化和物理實(shí)驗(yàn)的對(duì)比結(jié)果

由圖4可以看出來(lái)，數(shù)值模擬和物理實(shí)驗(yàn)在Fr＜1 時(shí)均表現(xiàn)出了非線(xiàn)性，拖曳力隨速度的增加先增大后減小。當(dāng)Fr＞1 時(shí)，拖曳力隨著速度的增大而增大。這一現(xiàn)象在拖曳系數(shù)的變化規(guī)律中更加明顯，而且當(dāng)Fr＞1 后，拖曳系數(shù)趨近于1，和單層流拖曳系數(shù)的變化規(guī)律相似，見(jiàn)圖5中數(shù)值模擬和物理實(shí)驗(yàn)的對(duì)比結(jié)果。

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

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

對(duì)于密度線(xiàn)性變化的連續(xù)分層流體，當(dāng)h=35 cm 時(shí)，其內(nèi)波相速度C：

而本研究采用兩層流體線(xiàn)性假設(shè)下的內(nèi)波相速度C：

圖6 數(shù)值模擬中流體分層條件和物理實(shí)驗(yàn)的對(duì)比

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

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

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

3.2 冰脊拖曳力隨入水深度的變化雙層流中冰脊拖曳力、拖曳系數(shù)隨無(wú)量綱入水深度B 的變化規(guī)律見(jiàn)圖8所示。

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

由拖曳系數(shù)的定義可知，拖曳系數(shù)和冰脊運(yùn)動(dòng)速度的平方呈反比。因此當(dāng)速度較小的時(shí)候，拖曳系數(shù)對(duì)拖曳力的變化更加敏感。拖曳力的非線(xiàn)性變化出現(xiàn)在弗勞德數(shù)較小的區(qū)間，此時(shí)拖曳系數(shù)隨Fr 的變化趨勢(shì)更為顯著（圖7）。因此圖8中Fr 較小時(shí)拖曳系數(shù)的變化幅度要大于Fr 較大時(shí)的變化幅度。

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

圖8 冰脊拖曳系數(shù)隨無(wú)量綱入水深度的變化結(jié)果

圖9 冰脊入水深度T=6 cm，速度V=15 cm/s 數(shù)值模擬的內(nèi)波

圖10 數(shù)值模擬內(nèi)波無(wú)量綱化波高和物理實(shí)驗(yàn)對(duì)比結(jié)果

圖11 冰脊入水深度T=8 cm，速度V=7 cm/s 數(shù)值模擬內(nèi)波形態(tài)和物理實(shí)驗(yàn)的對(duì)比結(jié)果（黑色虛線(xiàn)為初始分界面）

隨著弗勞德數(shù)的增長(zhǎng)，冰脊上游的擾動(dòng)消失，僅冰脊下方出現(xiàn)波谷狀態(tài)。此時(shí)，下游內(nèi)波仍處于比較穩(wěn)定的狀態(tài)。而當(dāng)Fr＞1 時(shí)，由于冰脊的拖曳速度超過(guò)了內(nèi)波的相速度，內(nèi)波擾動(dòng)僅向冰脊下游傳播，而且第一個(gè)波峰逐漸變得平坦，并逐漸遠(yuǎn)離冰脊。當(dāng)Fr=2.1（V=30 cm/s）時(shí)，內(nèi)波形態(tài)變得更加雜亂，波峰不斷遠(yuǎn)離冰脊，最終擾動(dòng)逐漸消失，流態(tài)過(guò)渡到超臨界區(qū)域。

其他入水深度條件下產(chǎn)生的內(nèi)波形態(tài)和T=6 cm 的情形一樣，但波峰高度有所不同。為進(jìn)行數(shù)值模擬和物理實(shí)驗(yàn)的比較，把內(nèi)波最低點(diǎn)與最高點(diǎn)的高度差定義為內(nèi)波波高H，H/h1為無(wú)量綱化波高，圖10是不同入水深度情況下無(wú)量綱波高的對(duì)比。

圖10中數(shù)值模擬的波高比物理實(shí)驗(yàn)結(jié)果偏大，但趨勢(shì)一致。在跨臨界狀態(tài)，波高隨著速度的增大而增大；在超過(guò)臨界狀態(tài)后，波高開(kāi)始趨于平穩(wěn)。數(shù)值模擬與物理實(shí)驗(yàn)的差異還是由于上文中提到的密度分層的不同，它導(dǎo)致了圖11中物理實(shí)驗(yàn)的摻混現(xiàn)象比較明顯，消耗了更多的能量。而數(shù)值模擬中的波形則比較穩(wěn)定，因此內(nèi)波波高相對(duì)較大。

3.4 內(nèi)波對(duì)冰脊拖曳力的影響由圖7中的對(duì)比可知，在拖曳力變化的非線(xiàn)性區(qū)間，內(nèi)波對(duì)冰脊拖曳力產(chǎn)生了重要影響，下面從垂直于冰脊運(yùn)行方向和平行于冰脊運(yùn)行方向兩個(gè)方面對(duì)該問(wèn)題進(jìn)行分析。首先對(duì)于垂直方向，統(tǒng)計(jì)了非線(xiàn)性區(qū)間的流體界面垂直偏移的最大值D，該距離隨弗勞德數(shù)的變化規(guī)律見(jiàn)圖12。從圖中可知，偏移距離D 和拖曳力F 在跨臨界區(qū)域相關(guān)性較高，相關(guān)系數(shù)R2＞0.94。分界面位移的變化可以反映出流體勢(shì)能的變化，同時(shí)，由于上層流體通道束窄，流速也會(huì)增加。同一流線(xiàn)在該位置勢(shì)能增加，動(dòng)能增加，而且下層流體對(duì)上層流體的剪切力做負(fù)功，流線(xiàn)的壓強(qiáng)勢(shì)能則會(huì)降低，在分界面位移最高處形成一個(gè)相對(duì)低壓區(qū)。

圖12 數(shù)值模擬分界面位移和冰脊拖曳力關(guān)聯(lián)圖（藍(lán)線(xiàn)是拖曳力，黑線(xiàn)是位移）

對(duì)于平行冰脊運(yùn)動(dòng)方向，冰脊拖曳力主要是由冰脊前后壓強(qiáng)場(chǎng)的變化引起的。因此對(duì)冰脊附近的壓強(qiáng)場(chǎng)進(jìn)行后處理，將流場(chǎng)總壓強(qiáng)減去靜壓，得到流場(chǎng)的動(dòng)壓強(qiáng)。這樣更容易反映出冰脊前后壓強(qiáng)的變化，結(jié)果見(jiàn)圖13。

圖13是不同運(yùn)動(dòng)速度冰脊附近的壓強(qiáng)場(chǎng)和流場(chǎng)。可以看到雙層流體中不僅會(huì)在冰脊后形成一個(gè)尾流漩渦場(chǎng)（和單層相同），而且還會(huì)在內(nèi)波擾動(dòng)的波峰、波谷處形成漩渦場(chǎng)。每一個(gè)漩渦場(chǎng)都是一個(gè)相對(duì)低壓區(qū)，渦動(dòng)中心的壓強(qiáng)最低，這也是冰脊產(chǎn)生拖曳力的原因。在單層流中，冰脊拖曳力只是受尾流漩渦場(chǎng)一個(gè)低壓區(qū)的影響，而且隨著速度的增大，漩渦中心的相對(duì)壓強(qiáng)降低。然而在雙層流中，冰脊拖曳力要受尾流漩渦場(chǎng)和冰脊后第一個(gè)內(nèi)波波峰漩渦場(chǎng)的雙重影響。尤其當(dāng)拖曳速度小于內(nèi)波相速度時(shí)，內(nèi)波擾動(dòng)向冰脊上下游同時(shí)傳播，此時(shí)內(nèi)波的形態(tài)相對(duì)冰脊比較穩(wěn)定，跟隨冰脊的運(yùn)動(dòng)而運(yùn)動(dòng)。波峰的漩渦場(chǎng)和尾流漩渦場(chǎng)共同作用，對(duì)冰脊拖曳力產(chǎn)生影響。

圖13 數(shù)值模擬分層流體冰脊附近動(dòng)壓強(qiáng)分布云圖（黑色實(shí)線(xiàn)為內(nèi)波波面）

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

4 冰脊底角的影響

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

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

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

圖15 T=6 cm，拖曳系數(shù)隨角度的變化：Fr≤0.71 為非線(xiàn)性區(qū)間、Fr＞0.71 為單調(diào)區(qū)間

5 結(jié)論

雙層流環(huán)境下冰脊拖曳力隨著冰脊運(yùn)動(dòng)速度和冰脊入水深度的變化而變化，對(duì)應(yīng)的冰脊拖曳系數(shù)受弗勞德數(shù)和冰脊的無(wú)量綱入水深度的影響。數(shù)值模擬得到的冰脊拖曳力和拖曳系數(shù)與物理實(shí)驗(yàn)結(jié)果基本吻合，在跨臨界區(qū)域內(nèi)都表現(xiàn)出冰脊拖曳力和拖曳系數(shù)的非線(xiàn)性變化規(guī)律。數(shù)值模擬產(chǎn)生的內(nèi)波形狀和物理實(shí)驗(yàn)結(jié)果一致，但由于分界面密度分布的差異，物理實(shí)驗(yàn)中存在更多的摻混現(xiàn)象，消耗了更多的內(nèi)波能量。

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

單個(gè)冰脊的局地拖曳系數(shù)對(duì)冰-水總拖曳系數(shù)的確定具有重要作用，結(jié)合流體分層的物理實(shí)驗(yàn)研究，本文利用數(shù)值模擬方法探討了雙層流中冰脊拖曳系數(shù)隨冰脊入水深度和弗勞德數(shù)變化關(guān)系，為分層流體冰-水拖曳系數(shù)的參數(shù)化奠定了基礎(chǔ)。除了冰脊入水深度和運(yùn)動(dòng)速度，冰脊形態(tài)和分層水深的變化也會(huì)對(duì)冰脊局地拖曳系數(shù)產(chǎn)生影響，需要進(jìn)一步的試驗(yàn)研究和理論分析，從而發(fā)展完整的冰脊拖曳系數(shù)參數(shù)化方案。