基于Rankine源法的氣墊船破冰數值模擬

李宇辰, 劉巨斌, 丁志勇, 張志宏

(海軍工程大學 理學院，武漢 430033)

基于Rankine源法的氣墊船破冰數值模擬

李宇辰, 劉巨斌, 丁志勇, 張志宏

(海軍工程大學 理學院，武漢 430033)

基于黏彈性薄板假設和勢流理論，建立了氣墊船在水面、冰面、碎冰面上航行時引起冰層-水層振動問題的統一的理論數學模型，提出了采用Rankine源與有限差分相結合的數值計算方法，通過C語言與Matlab聯合編程，建立了冰層-水層的位移響應、應力分布以及破裂效果的預報方法。針對黃河水域冰層特點，對氣墊船以亞臨界、臨界以及超臨界速度航行時冰層的應力分布以及破裂情況進行了數值模擬。結果表明：所建立的理論模型和計算方法可以反映冰層的位移變形及其內部應力的分布特征，可以捕捉到氣墊船破冰的臨界速度，計算得到的冰層振動特性與已有文獻的研究結果一致。

氣墊船；破冰；Rankine源法；有限差分法

我國黃河寧蒙段地理環境特殊，每年春季開河時，水流因冰層阻礙常在河段的彎曲和狹窄處引發河水上漲，嚴重時將造成凌汛災害，致使人民生命財產遭受重大損失，如何及時破冰除險成為每年防災減災的重要任務。黃河寧蒙段河道狹窄、淺灘迭出，難于利用破冰船實施破冰作業，目前主要以爆破方式為主。防凌實踐經驗表明，采用爆破方式破冰不僅需要疏散當地群眾、耗費大量人力物力，而且還存在破壞河道生態環境以及安全隱患等問題，因此亟需開發安全、環保、可靠、高效的破冰新方法。

氣墊船在冰層上以臨界速度航行時可以激勵冰層大幅變形而達到破冰的效果。Squire等[1]對移動載荷在冰層上的運動理論以及實驗進行了總結。近年來，俄羅斯對氣墊船破冰的原理進行了研究[2-6]。之后國內也開展了相關研究工作：文獻[7-8]采用積分變換方法求解了淺水岸壁條件下脈沖載荷和移動載荷激勵冰層的位移響應問題，并導出了淺水臨界速度的近似計算公式；文獻[9]采用LS-DYNA動力學分析軟件對氣墊船破冰過程進行了數值模擬，并根據計算結果分析了氣墊船的破冰機理。上述研究工作主要針對簡單邊界條件下氣墊載荷激勵冰層變形的響應問題，對氣墊載荷在水面、冰面、碎冰面上航行引起的冰層振動、冰層內應力分布以及冰層破裂情況等問題研究很少。本文在文獻[10-11]基礎上，采用Rankine源與有限差分相結合的數值計算方法，進一步對氣墊船航行引起的冰層變形、應力分布以及破裂情況進行數值模擬，為氣墊船在黃河實施有效破冰提供技術手段。

1 理論模型

為建立統一的冰-水振動模型，設半無限大浮冰層覆蓋于水面之上(冰厚為h，冰密度為ρ1)，另一半為自由表面(水深為H，水密度為ρ2)，自由表面與冰面的分割線是直線，氣墊船沿冰-水分割線方向在冰面或水面上以勻速U作直線運動。建立隨船運動的坐標系oxyz，oxy與未擾動水面重合，x方向與氣墊船運動方向一致，z軸垂直向上，z=0為未擾動水面，z=-H為水底，如圖1所示。

(a)正視圖

(b)俯視圖圖1 理論模型示意圖Fig.1 Sketch of the theoretical model

假設冰層為各向同性、勻質、厚度均勻的黏彈性薄板，水為理想不可壓縮流體作無旋運動，速度勢為Φ(它可以分解為來流速度勢-Ux與擾動速度勢φ(x,y,z)之和)，則冰層或自由表面垂向位移的控制方程可寫為

D(1-τU?/?x)▽4w+pA-pW+

ρ1gh+ρ1hU2?2w/?x2=0

(1)

式中：w為冰層或水面的垂向位移；D=Eh3/[12(1-μ2)]為冰層彎曲剛度，其中E為彈性模量，μ為泊松比；g為重力加速度；τ為冰層的松弛時間；pW為冰-水交界面處水的壓強；pA為自由表面或冰面上氣體的壓強； ▽4=(▽2)2=(?2/?x2+?2/?y2)2。

擾動速度勢在水域內滿足拉普拉斯方程：

?2φ/?x2+?2φ/?y2+?2φ/?z2=0

(2)

在水底z=-H處滿足不可穿透條件：

?φ/?z=0

(3)

在自由表面上滿足動力學邊界條件：

ρ2gw-ρ2U?φ/?x=-pAe

(4)

式中：pAe=pA-pA∞為自由表面的相對壓強，其中pA∞為無窮遠處自由表面上的壓強。

在冰-水交界面z=0處滿足運動學條件：

?φ/?z=-U?w/?x

(5)

在冰-水交界面z=0處滿足動力學邊界條件：

pW=pA∞+ρ1gh-ρ2gw+ρ2U?φ/?x

(6)

無窮遠處邊界條件為

w=0,▽w=0

(7)

將式(5)～式(7)與式(1)進行耦合可得：

(8)

至此，式(8) 和式(2)～式(4)一起構成了求解氣墊船激勵冰-水振動的統一數學模型。

當h=0時，式(8)退化為式(4)，則式(2)～式(4)可用于求解純水面興波問題。

當h≠0，D=0時，式(8)退化為

(9)

此時式(2)、式(3)以及式(9)可用于求解碎冰問題。

將氣墊壓強分布寫作雙曲正切函數的乘積形式，即：

(10)

式中：p0是氣墊船特征壓強；L是氣墊船的長度；B是氣墊船的寬度；(x0,y0)是氣墊船的中心位置；α和β是反映氣墊壓強變化程度的參數。

在求得擾動速度勢φ后，水域自由表面波高或冰層垂向位移可由式(11)進行計算：

(11)

氣墊船興波阻力的計算公式為

(12)

式中：S為氣墊壓強作用區域。

冰層內的法向應力σxx、σyy和切向應力σxy為

(13)

綜上所述，求解擾動速度勢φ是求解冰層變形、冰層內應力分布以及氣墊船航行興波阻力的關鍵所在。下面對擾動速度勢φ的計算方法進行說明。

2 數值計算方法

本文采用Rankine源法與有限差分法相結合的數值計算方法來求解擾動速度勢φ。具體做法是，在自由表面上方一定距離、一定區域的水平面上分布面元基本解(Rankine源)，因此拉普拉斯方程式(2)自動滿足，水底邊界條件式(3)通過鏡像法得到滿足，而控制方程式(8)則用來確定源強。

2.1 面元法

將擾動速度勢表示為面元基本解的線性組合：

(14)

圖2繪制了基本解和控制點的布置方式。在自由表面上方一定高度處的水平面上分布n=nx×ny個矩形面元，矩形的中心在水面上的投影為控制點，因而控制點的個數等于面元的個數。

(a) 正視圖

(b) 俯視圖圖2 面元和控制點的布置Fig.2 Distribution of panels and control points

2.2 控制方程的離散方法

由于采用面元法以及鏡像法使得該問題的邊界條件以及拉普拉斯方程自動得到滿足，因此求解水波-冰層變形問題的關鍵在于處理式(8)。該式中除?/?z使用解析方法運算外，其余偏導數運算均采用數值微分的方法進行計算，其中?/?x方向采用迎風差分。具體離散格式參見文獻[10]。

通過處理式(8)并聯立式(14)可以得到關于σj的線性方程組，求解該方程組可以得到源強分布，進而可以得到擾動速度勢φ。

3 算法實現與結果分析

3.1 算法實現

基于理論模型與數值計算方法，采用C語言與Matlab聯合編程的方法編寫數值計算源程序。程序中源強的確定部分采用OPENMP實現并行計算，提高了程序的運行效率。

3.2 計算結果與分析

數學模型與計算方法的正確性已經過有效驗證，驗證結果表明：純水面、純碎冰面以及純冰面情況下，興波阻力系數的數值計算結果與理論計算值符合良好；不同速度下冰層的最大下陷位移計算結果不僅與Takizawa等[13]的實驗結果相一致，且能捕捉到氣墊船的第一臨界速度和第二臨界速度。下面參照我國黃河實際冰層情況(見表1)，以純冰面情況為例給出冰層變形、應力分布以及冰面破裂計算結果。

表1 冰層物理參數Tab.1 Physical parameters of ice sheet

假設氣墊船長L=20 m，寬B=10 m，水深H=3 m，冰厚h=0.2 m，τ=0.5，氣墊壓力p0=3 000 Pa，計算域為16L×16L，氣墊船中心位于(9.6L,8L)處。根據表1冰層的物理參數，采用張志宏等給出的淺水臨界速度近似公式計算得到臨界速度約為5.5 m/s，以此為參考，對氣墊船以亞臨界速度(4 m/s)、臨界速度(5.5 m/s)和超臨界速度(10 m/s)航行時引起的冰層變形、冰層應力分布以及冰層破裂情況進行數值計算，結果如圖3，4所示。

圖3為不同速度下冰層的位移變形，圖4為相應速度下的應力等值線圖(以σxx為例)。從圖中可以看出：當氣墊船速度為4 m/s(亞臨界速度)時，冰層位移變形很小，氣墊船前后波形幾乎呈對稱分布，此時與靜態載荷作用下的冰層變形情況類似；當氣墊船速度為5.5 m/s(臨界速度)時，氣墊船前方出現橫波，主要為周期較短幅值較小衰減較快的彈性波，后方為橫波和散波，主要為周期較長幅值較大衰減較慢的重力波，冰層最大下陷位移幅值相比于亞臨界速度情況有大幅增加，氣墊船船首接近于波峰，船尾接近于波谷；當氣墊船速度進一步增大為10 m/s(超臨界速度)時，氣墊船前方仍然存在彈性波，而其后方為“V”型重力波，此時冰層最大下陷位移和應力分布相比于臨界速度時要大幅減小。

(a) 4 m/s(亞臨界速度)

(b) 5.5 m/s(臨界速度)

(c) 10 m/s(超臨界速度)圖3 不同速度時的冰層位移響應Fig.3 displacement response of ice sheet at different speeds

(a) 4 m/s(亞臨界速度)

(b) 5.5 m/s(臨界速度)

(c) 10 m/s(超臨界速度)圖4 冰層應力等值線圖(σxx)Fig.4 Contour of stress on ice sheet (σxx)

如果以5×105Pa作為冰層的破壞應力，那么將超過該臨界值的面元標記為黑色則可以得到如圖5所示的冰層破裂效果。由圖可見，在氣墊船航速較小時，冰面除氣墊船附近以外極少有破裂面元；而當氣墊船以臨界速度航行時，冰面上出現大范圍的破裂面元；當氣墊船以超臨界速度航行時，冰面上的破裂單元相對于臨界速度時又大幅減少。從上述變化規律可以看出，該數值計算程序可以捕捉到氣墊船航行的臨界速度以及不同航速條件下的破冰效果。氣墊船以臨界速度航行時破冰效果最好，如圖5(b)所示。

(a) 4 m/s(亞臨界速度)

(b) 5.5 m/s(臨界速度)

(c) 10 m/s(超臨界速度)圖5 氣墊船破冰示意圖Fig.5 Sketch of ice-breaking with ACV

氣墊船破冰的機理是移動載荷以臨界速度航行時，移動載荷作用于冰層-水層系統的波動能量可以不斷累積，激勵冰層引起聚能共振增幅效應，從而誘發冰層大幅變形，并導致冰層內部的拉壓、彎曲應力超過其屈服極限而破裂。破冰效果預測與所采用的破冰準則有關，若采用文獻[12]提出的關系式作為破冰準則，即：λ=max|?w/?x|≥0.04，則可以得到如圖6所示的破冰效果。此時可以看出，在氣墊船以亞、超臨界速度航行時破冰效果差，而以臨界速度航行時仍能得到與圖5(b)一樣好的破冰效果。

(a) 4 m/s(亞臨界速度)

(b) 5.5 m/s(臨界速度)

(c) 10 m/s(超臨界速度)圖6 氣墊船破冰示意圖Fig.6 Sketch of ice-breaking with ACV

4 結 論

對氣墊船激勵冰層-水層變形問題建立了統一的理論數學模型和數值計算方法。針對黃河水域冰層特點，對氣墊船在純冰面上航行引起的冰層變形、應力分布以及破裂情況進行了數值計算。計算結果可以反映冰層的位移變形及其內部應力的分布特征，可以捕捉到引起冰層大幅變形的移動載荷臨界速度，并可根據給定的破冰準則預測氣墊船在不同航速下的破冰效果。本文研究工作可以推廣用于更復雜工況下(如水、冰以及碎冰兩兩任意組合)移動氣墊載荷激勵冰層-水層響應問題的數值預報。破冰準則的合理選擇還需要進一步優化和通過實驗的驗證。

[1] SQUIRE V A, HOSKING R J, KERR A D, et al. Moving loads on ice plates[M]. The Netherlands: Kluwer Academic Publishers, 1996.

[2] KOZIN V M, MILOVANOVA A V. The wave resistance of amphibian aircushion vehicles in broken ice[J]. Journal of Applied Mechanics and Technical Physics, 1996,37(5):634-637.

[3] KOZIN V M, POGORELOVA A V. Effect of broken ice on the wave resistance of an amphibian air-cusion vehicle in nonstationary motion[J]. Journal of Applied Mechanics and Technical Physics, 1999,40(6):1036-1041.

[4] KOZIN V M, POGORELOVA A V. Variation in the wave resistance of an amphibian air-cusion vehicle moving over a broken-ice land[J]. Journal of Applied Mechanics and Technical Physics,2007，48(1):80-84.

[5] KOZIN V M, POGORELOVA A V. Wave resistance of amphibian aircushion vehicles during motion on ice fields[J]. Journal of Applied Mechanics and Technical Physics, 2003,44(2):193-197.

[6] ZHESTKAYA V D,KOZIN V M. Stress-deformed state of a semi-infinite Ice sheet under the action of a moving load[J]. Journal of Applied Mechanics and Technical Physics, 1994,35(5):745-749.

[7] 鹿飛飛, 張志宏, 胡明勇,等. 淺水岸壁條件下脈沖載荷引起的黏彈性浮冰層位移響應[J]. 振動與沖擊, 2015, 34(14)：142-146.

LU Feifei, ZHANG Zhihong, HU Mingyong, et al. Displacement response of viscoelastic floating ice sheet subjected to impulse load under different bank conditions[J]. Journal of Vibration and Shock, 2015, 34(14)：142-146.

[8] 張志宏, 鹿飛飛, 丁志勇,等. 勻速移動載荷激勵浮冰層大幅響應的臨界速度[J]. 華中科技大學學報(自然科學版), 2016，44(2):107-111.

ZHANG Zhihong, LU Feifei, DING Zhiyong, et al. Critical speed of a sharp response for floating ice sheet subjected to moving load with uniform speed[J]. Journal of Huazhong University of Science and Technology(Natural Science Edition), 2016, 44(2):107-111.

[9] 盧再華, 張志宏, 胡明勇,等. 全墊升式氣墊船破冰過程的數值模擬[J]. 振動與沖擊, 2012, 31(24):148-154.

LU Zaihua, ZHANG Zhihong, HU Mingyong, et al. Numerical simulation for ice-breaking process of an amphibian air cushion vehicle[J]. Journal of Vibration and Shock, 2012, 31(24):148-154.

[10] 劉巨斌, 張志宏, 張遼遠, 等. 氣墊船興波破冰問題的數值計算[J]. 華中科技大學學報(自然科學版), 2012, 40(4):91-95.

LIU Jubin, ZHANG Zhihong, ZHANG Liaoyuan, et al. Numerical computation of broken ice by air-cushion vehicles in wave making[J]. Journal of Huazhong University of Science and Technology(Natural Science Edition),2012,40(4):91-95.

[11] 劉巨斌, 張志宏, 張遼遠, 等. 邊界元—有限差分法在氣墊船破冰數值模擬中的應用[J]. 海軍工程大學學報, 2013, 25(3):50-55.

LIU Jubin, ZHANG Zhihong, ZHANG Liaoyuan, et al. Application of mixed BEM and FDM in numerical simulation of ice-breaking by air cushion vehicle[J]. Journal of Naval University of Engineering, 2013, 25(3):50-55.

[12] POGORELOVA A V, KOZIN V M. Motion of a load over a floating sheet in a variable-depth pool[J]. Journal of Applied Mechanics and Technical Physics, 2014,55(2):335-344.

[13] TAKIZAWA T. Deflection of a floating sea ice sheet induced by a moving load[J]. Cold Regions Science and Technology, 1985,11:171-180.

NumericalsimulationofACVs’ice-breakingbasedonRankinesourcemethod

LI Yuchen, LIU Jubin, DING Zhiyong, ZHANG Zhihong

(College of Science，Naval University of Engineering，Wuhan 430033，China)

Based on the viscoelastic thin plate assumption and the potential flow theory, a unified theoretical dynamic model was established to study ice layer-water layer vibration problems due to ACV sailing on water surface, ice surface or broken ice surface. The numerical method combing Rankine source method and the finite difference method (FDM) was proposed, the prediction method for ice layer-water layer displacement responses, stress distribution and ice-breaking effect was established using the simulation program written with C language and Matlab. Aiming at characteristics of ice layer in Yellow River water field, its stress distribution and ice-breaking effect were simulated numerically when ACVs sailing at subcritical speed, critical one and supercritical one. These simulated results showed that the established theoretical model and simulation method can not only reflect ice layer displacement deformation and its internal stress distribution characteristics, but also capture the critical speed of ACVs’ ice-breaking; the ice layer’s vibration features calculated agree well with those in the existing literature.

air cushion vehicle (ACV); ice-breaking; Rankine source method; FDM

國家自然科學基金(51479202)

2016-08-08 修改稿收到日期：2016-09-07

李宇辰 男,博士生,1988年生

張志宏 男,博士,教授,博士生導師,1964年生

O35

A

10.13465/j.cnki.jvs.2017.23.005