朱仁慶, 張曦, 李志富 (江蘇科技大學 船舶與海洋工程學院,江蘇 鎮江 212003)
隨著全球氣候變暖,海冰冰厚變薄使得北極資源的開發以及航線開通成為可能。許多國家競相儲備開發北極資源和北極航道,準備隨時進駐北極。因此,合理認識和利用冰緣區內波-冰作用規律,有利于加強對冰區內波浪傳播特性的了解和提高極地資源利用率具有重要意義。對于浮冰與波浪的相互作用,Harms[1-2]通過試驗回歸出恒密度厚度的二維浮冰漂移運動估算公式。當浮冰的橫搖周期略大于波浪周期時,浮冰具有最大漂移速度。Lever等[3]在水池試驗中研究了冰山自身形狀對其在波浪中運動的影響。Huang等[4]研究發現浮冰漂移速度隨著橫搖周期與波浪周期的比值增大而增大。Bennetts等[5]在試驗中發現浮冰的存在阻礙了波浪傳播,且隨著波陡的增大,浮冰對波浪的阻礙作用增強。Yiew等[6]發現當波陡較大時,波浪沖洗現象對浮冰的運動響應有明顯抑制作用。郭春雨等[7]試驗研究了浮冰在波浪中的縱向運動,同時發現了波浪沖洗現象。波浪在冰區的傳播比例隨著波長和波頻率的減小而增大,且波浪傳播比例增大也增強了其對浮冰的破壞能力[8]。
由于模型試驗具有成本過高、存在設備干擾的缺點,數值模擬方法被越來越多的用于研究波-冰作用。其中,勢流分析方法因其成本低、計算速度快等優點受到了學者的青睞。Shen等[9]通過理論研究了反射波、阻力系數和附加質量影響浮體漂移運動的的因素,如波的反射、附加質量等。Meylan[10]研究發現浮冰運動很大程度上依賴于浮冰形狀,波浪散射主要取決于入射波方向。Bennetts等[11]研究了非均質浮冰對波浪的散射。Smith等[12]對均質浮冰對波浪散射現象求解進行了研究,并在此方法的基礎上發展了非均質浮冰對波浪散射的求解。Meylan等[13]通過改進莫里森方程研究大波長規則波中浮冰的縱蕩運動,與試驗結果具有較好的一致性。Toffoli等[14]提出了波浪海冰相互作用的試驗模型來驗證冰緣區波浪衰減的理論模型。由于波-冰作用中固有的強非線性現象無法通過勢流理論精確計算出,因此計算流體力學動態模擬技術的迅速發展研究波-冰作用起到了重要的作用。Bai等[15]研究表明粘流軟件比勢流軟件更適合于模擬浮冰運動。Huang等[16]通過粘流軟件模擬了波浪沖洗冰體表面和波浪散射現象。
在浮冰與波浪的相互作用中,浮冰在波浪作用下會發生非形變運動。與此同時,波浪場也會由于浮冰的運動而發生改變。目前關于波-冰作用的研究對于相互影響已經有了充分了解,但對波浪沖洗現象對浮冰運動的影響研究較少。波浪沖洗冰體表面現象雖有分析,但是很少對波浪沖洗現象、浮冰運動以及波浪參數之間定性定量關系的研究。本文通過粘流模型數值模擬了波浪中浮冰的運動,利用頻譜分析研究了浮冰運動響應。分析了波浪沖洗現象對浮冰運動的影響,對比在有/無波浪沖洗現象工況下的浮冰運動。并研究了不同波浪參數下波浪沖洗現象對浮冰垂蕩、縱搖運動的影響。
浮冰區域流體運動滿足連續方程和動量方程。考慮到湍流流動,通常采用時均法將流體運動控制方程中的各物理量分解為2部分,分別為時間平均值和相對于平均值的脈動值,得到雷諾時均方程組,其表達式為:

(1)
(2)

Pk-ρ(ε-ε0)+Sk
(3)
(4)

本文通過速度入口造波方式生成波浪,通過在入口處和出口處分別通過松弛區域消波技術和阻尼消波方法實現消波從而減少壁面反射。本文選取斯托克斯五階波浪[18-19],使模擬的波浪更符合實際的規則波。在造波邊界處的速度和波面瞬時升高滿足以下條件:
x方向速度:
(5)
z方向速度:
(6)
波面瞬時升高:
(7)
其中各系數分別為:
式中:ω、d、k分別是圓頻率、水深和波數;c=cosh(kd);λ為常系數,其余各項系數參考文獻[20]。
入口消波區的消波原理是通過在動量方程中添加源項qφ強迫波浪在速度入口附近轉換成入射波:
qφ=-γρ(φ-φ*)
(8)
式中:γ為強迫系數;φ為數值模擬中的動量輸運方程的結果;φ*為原動量方程的理論解。冰水相互作用的數值水池模型如圖1所示。

圖1 冰水相互作用數值水池Fig.1 Diagram of the numerical wave tank
出口消波區的原理是通過在Z軸方向上添加阻尼完成消波作用,其中阻尼項為:
(9)
式中:xsd為起始消波位置;xed為結束消波位置;f1、f2、nd為阻尼消波的參數;w為Z方向的速度。
在隨浮冰平移和旋轉的運動坐標系中,浮冰的運動方程為:
(10)
(11)
式中:m是浮冰質量;V是浮冰運動速度;Ω是浮冰旋轉角速度;I是浮冰轉動慣量;M是浮冰受到的合外力矩;f是浮冰受到的的合外力。
2.1.1 波浪驗證
根據極地冰緣區出現的波浪高度和頻率分布規律和波浪參數與浮冰特征尺寸的比例關系,在相同的波浪高度下,對浮冰在一系列波長(波陡偏大)波浪中的運動進行了數字模擬和分析研究,工況1~5波高恒定為0.04 m,并且入射波長為0.6、0.7、0.8、0.9、1.0 m。
液面波高探測點位置選取為圓柱所在位置,其中選取工況5的波浪來記錄探測點處液面升高過程,將斯托克斯五階波理論解與其監測到的波高時歷曲線數據進行比較。對于波陡小的波浪,造波效果與其設置的波高H方向的網格尺寸Δz有重要的關聯,而能否生成目標波浪決定了計算結果的準確與否。因此需要再進行收斂性分析,這里分別選取3組不同大小的自由液面處波高方向網格和3組不同時間步進行生成波浪進行分析,Δz分別為H/5、H/10及H/15。波長方向的網格尺寸Δx是波高方向網格尺寸Δz的4倍。網格與時間步設置如表1所示。

表1 不同網格尺寸與時間步設置Table 1 Different mesh size and time step settings
如圖2所示,數值水池生成的波浪在頻率、波高等波浪參數在網格B上的誤差最小,與斯托克斯五階波理論解吻合較好。因此,波高方向上網格尺寸Δz為H/10被選取進行數值模擬。將該套網格在3種不同時間步下進行了計算,如圖3所示,數值模擬時間步Δt為0.002 s的工況的結果收斂。因此,出于效率和資源上的綜合考慮,波高方向網格尺寸Δz為H/10,時間步長Δt為0.002 s被選入進行后面的數值模擬。

圖2 波高方向不同網格尺寸Fig.2 Dimensions of different mesh numbers in the wave height direction

圖3 不同時間步工況與理論值對比Fig.3 Comparison of wave conditions and theoretical solution values at different time steps
2.1.2 模型驗證
本文選取文獻[6]中單個浮冰與規則波相互作用試驗研究,提取試驗中浮冰的垂蕩幅值時歷曲線來驗證數值模擬計算結果。浮冰模型為圓柱體,半徑為0.2 m,厚度為0.015 m,密度為636 kg/m3,浮冰邊緣附有薄壁結構,薄壁結構高為0.05 m,厚度為0.025 m,在浮冰周圍施加彈簧力以防止浮冰產生大幅度漂移。在保證無池壁反射的前提下,適當縮減了水池的長度和寬度,數值模型中水池寬2 m,長7 m,水深0.83 m。入射波頻率為1.25 Hz,波幅為0.085 m。浮冰中心距離造波入口4 m處。計算模型中網格劃分遵循上述網格參數設置,在浮冰周圍加密以保證計算精度,網格劃分圖如圖4所示,其中對于k-ε模型需要滿足30≤y+≤300,一般以接近30為佳,因此y+取值為30,第1層網格厚度為0.000 47 m。

圖4 重疊網格劃分Fig.4 The arrangement of the mesh
由圖5可知,試驗值與數值模擬中的浮冰垂蕩運動基本一致,浮冰垂蕩運動峰值以及到達相鄰峰值所需的周期大致吻合,由于偏差較小,數值模擬結果相較于實驗結果沒有明顯偏高或偏低趨勢,因此可以將偏差視為數值計算或實驗測量精度的判斷依據,說明本文采用的數值模擬方法是合理和準確的。

圖5 浮冰垂蕩運動對比Fig.5 The comparison of ice floe heave motion
如圖6所示,數值水池尺寸為4.5 m×2 m×0.8 m,水深0.5 m。數值水池的前后邊界設置為對稱邊界條件;左側邊界、上側邊界和下側邊界設置為速度入口邊界條件;右側邊界設置為壓力出口邊界條件。入口以及出口消波區區域大小均取為一個波長。

圖6 波-冰作用的三維數值水池Fig.6 Three-dimensional numerical tank of wave-ice interaction
將浮冰視為剛體,由于實際海況上浮冰互相之間經常會受到碰撞和波浪的沖洗,浮冰邊角的磨損使其形狀近似于圓柱,如圖7所示。將浮冰模型進行處理便于計算,故簡化為圓柱體和帶一定高度薄壁的圓柱體2種形狀,浮冰中心與入口處的距離是1.5 m。其中圓柱參數:直徑0.5 m;厚度0.05 m;吃水0.03 m;密度600 kg/m3。為對無波浪沖洗現象進行研究,將圓柱浮冰周圍增設高為0.05 m薄壁,研究表明薄壁對浮冰運動的影響可忽略不計[21]。

圖7 2種浮冰模型Fig.7 The geometry of two kind of ice floes
本文通過重疊網格技術模擬浮冰運動,在浮體運動區域以及自由液面區域進行網格加密,如圖8為重疊網格。

圖8 重疊網格Fig.8 Overset mesh
波浪的參數對浮冰運動響應有著關鍵的影響。為了進一步研究浮冰與入射波的相互作用,通過頻譜分析研究了不同波長下浮冰垂蕩η3、縱搖η5運動響應。
浮冰在波浪作用下會發生大幅的搖蕩運動,海水沖洗浮冰表面的現象頻頻出現。由數值模擬可以看出,波長增加,波浪沖洗現象會隨之減弱。圖9為在不同波長條件下一個周期的浮冰表面的上浪現象。圖9(a)是波長為0.6 m時波浪沖洗狀態,浮冰表面上浪嚴重;圖9(b)是波長為1.0 m時的波浪沖洗狀態,浮冰表面未完全被水體覆蓋,與圖9(a)相比,其所受的波浪沖洗現象較弱。

圖9 不同波長的波浪沖洗冰體表面過程Fig.9 Overwash process of waves with different wavelength
圖10為λ=0.6 m時的浮冰搖蕩運動頻譜分析圖。在相同的波浪激勵頻率作用下,沒有受到波浪沖洗的浮冰在運動幅值上要大于受到波浪沖洗的浮冰。產生上述差異的主要原因是在波浪與浮冰相互作用中存在波浪沖洗冰體現象,一方面堆積在冰體表面的水阻礙浮冰的運動;另一方面冰體表面的水體會產生碰撞繼而引起水體破碎,使波浪產生一定程度的衰減,使得浮冰接受的波浪能量較少。

圖10 波長λ=0.6 m的浮冰運動幅值對比Fig.10 Comparison of the motion amplitudes of the two ic floes, the wave condition was λ=0.6 m
圖11為不同波長下的波浪沖洗現象對浮冰的垂蕩和縱搖的影響。浮冰的垂蕩與縱搖隨著波長減小而減小。當波長較短時,波陡較大,波浪沖洗冰體表面現象較嚴重。由于堆積在冰體表面的水體較多以及波浪能量衰減較大,一定程度上阻礙了浮冰運動,使其運動幅度減弱。

圖11 無護欄浮冰在不同波長下的運動響應對比Fig.11 Comparison of the motion response of floating ice without guardrail at different wavelengths
由此可見,在給定波高情況下,波長越小,即波陡越大,波浪沖洗現象越嚴重。浮冰運動在波浪沖洗現象作用下受到了一定的抑制,波浪沖洗現象越強,其對浮冰的抑制作用越強。
通過模擬加薄壁的浮冰在波浪中的運動來研究波-冰作用。通過數值模擬可以看出,短波長下波-冰作用更加劇烈。浮冰的存在一定程度上阻礙了波浪的運動,波浪場存在較明顯的散射現象。
圖12(a)對應波長較短的工況,波浪場產生了明顯的波浪散射現象,浮冰前側和后側的波浪液面形狀由于散射的作用發生了明顯的擾動現象。圖12(b)對應波長較長的工況,波浪散射現象較弱,浮冰周圍波形沒有明顯變化。

圖12 有護欄浮冰在不同波長下與波浪相互作用Fig.12 Floes with guardrails interact with waves at different wavelengths
如圖13所示,在波浪激勵頻率下,圖13(a)中波長越小,浮冰的垂蕩運動幅值也越小。當波長較短時,浮冰大大阻礙了波浪的傳遞,同時也對波浪場造成了一定散射。由于波浪的散射消耗了一定的波能,使得浮冰的運動響應和受到的波浪力也隨之減弱。隨著波長增加,浮冰相對波長的尺度減小,進而浮冰對波浪的作用也減小,浮冰接收的波浪能量變大,因而波浪激勵頻率下的浮冰運動也隨之劇烈。圖13(b)中,浮冰的縱搖運動隨著波長的減小愈發劇烈。然而當波長變化至0.8 m時,波-冰特征長度比值為1.6,浮冰的縱搖運動開始減弱。這是因為波浪散射的程度隨著波陡的增加愈發嚴重,波能耗散也愈發嚴重,從而使浮冰受到的作用力和運動響應減弱。

圖13 有護欄浮冰在不同波長下的運動對比Fig.13 Comparison of the movement amplitude of ice floes with guardrails at different wavelengths
綜上所述,隨著波陡增大,波浪的散射現象越來越明顯。波浪散射現象對浮冰垂蕩具有一定程度的抑制作用。研究發現波浪散射現象越強,浮冰垂蕩受到的抑制作用也越強。其中,對于浮冰的縱搖運動,隨著波陡的增大,縱搖運動幅值增大,當波陡增大到一定值時,浮冰縱搖運動隨著波陡的減小而減小。
1)通過對比有沖洗/無沖洗浮冰的運動發現,波浪沖洗現象對浮冰的運動具有明顯的抑制作用。
2)隨著波陡變大,波浪沖洗現象和波浪散射現象增強,波浪沖洗現象對浮冰運動的抑制作用增強。同樣地,波浪散射現象對浮冰垂蕩運動具有抑制作用,波長越小,浮冰垂蕩運動受到的抑制作用越強。
3)波陡越大,浮冰的縱搖運動越大,當波陡增大到一定值時,由于波浪能量耗散嚴重,浮冰縱搖運動幅值隨波陡增大呈現出下降趨勢。
本文計算方法和分析結論對于進一步研究浮冰在波浪下的運動以及波浪中浮冰與結構物的相互作用具有參考價值。