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

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

基于Rankine源法的氣墊船破冰數(shù)值模擬

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

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

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

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

我國黃河寧蒙段地理環(huán)境特殊，每年春季開河時(shí)，水流因冰層阻礙常在河段的彎曲和狹窄處引發(fā)河水上漲，嚴(yán)重時(shí)將造成凌汛災(zāi)害，致使人民生命財(cái)產(chǎn)遭受重大損失，如何及時(shí)破冰除險(xiǎn)成為每年防災(zāi)減災(zāi)的重要任務(wù)。黃河寧蒙段河道狹窄、淺灘迭出，難于利用破冰船實(shí)施破冰作業(yè)，目前主要以爆破方式為主。防凌實(shí)踐經(jīng)驗(yàn)表明，采用爆破方式破冰不僅需要疏散當(dāng)?shù)厝罕姟⒑馁M(fèi)大量人力物力，而且還存在破壞河道生態(tài)環(huán)境以及安全隱患等問題，因此亟需開發(fā)安全、環(huán)保、可靠、高效的破冰新方法。

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

1 理論模型

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

(a)正視圖

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

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

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

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

(1)

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

擾動(dòng)速度勢在水域內(nèi)滿足拉普拉斯方程：

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

(2)

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

?φ/?z=0

(3)

在自由表面上滿足動(dòng)力學(xué)邊界條件：

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

(4)

式中：pAe=pA-pA∞為自由表面的相對壓強(qiáng)，其中pA∞為無窮遠(yuǎn)處自由表面上的壓強(qiáng)。

在冰-水交界面z=0處滿足運(yùn)動(dòng)學(xué)條件：

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

(5)

在冰-水交界面z=0處滿足動(dòng)力學(xué)邊界條件：

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

(6)

無窮遠(yuǎn)處邊界條件為

w=0,▽w=0

(7)

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

(8)

至此，式(8) 和式(2)～式(4)一起構(gòu)成了求解氣墊船激勵(lì)冰-水振動(dòng)的統(tǒng)一數(shù)學(xué)模型。

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

當(dāng)h≠0，D=0時(shí)，式(8)退化為

(9)

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

將氣墊壓強(qiáng)分布寫作雙曲正切函數(shù)的乘積形式，即：

(10)

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

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

(11)

氣墊船興波阻力的計(jì)算公式為

(12)

式中：S為氣墊壓強(qiáng)作用區(qū)域。

冰層內(nèi)的法向應(yīng)力σxx、σyy和切向應(yīng)力σxy為

(13)

綜上所述，求解擾動(dòng)速度勢φ是求解冰層變形、冰層內(nèi)應(yīng)力分布以及氣墊船航行興波阻力的關(guān)鍵所在。下面對擾動(dòng)速度勢φ的計(jì)算方法進(jìn)行說明。

2 數(shù)值計(jì)算方法

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

2.1 面元法

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

(14)

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

(a) 正視圖

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

2.2 控制方程的離散方法

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

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

3 算法實(shí)現(xiàn)與結(jié)果分析

3.1 算法實(shí)現(xiàn)

基于理論模型與數(shù)值計(jì)算方法，采用C語言與Matlab聯(lián)合編程的方法編寫數(shù)值計(jì)算源程序。程序中源強(qiáng)的確定部分采用OPENMP實(shí)現(xiàn)并行計(jì)算，提高了程序的運(yùn)行效率。

3.2 計(jì)算結(jié)果與分析

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

表1 冰層物理參數(shù)Tab.1 Physical parameters of ice sheet

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

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

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

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

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

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

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

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

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

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

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

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

氣墊船破冰的機(jī)理是移動(dòng)載荷以臨界速度航行時(shí)，移動(dòng)載荷作用于冰層-水層系統(tǒng)的波動(dòng)能量可以不斷累積，激勵(lì)冰層引起聚能共振增幅效應(yīng)，從而誘發(fā)冰層大幅變形，并導(dǎo)致冰層內(nèi)部的拉壓、彎曲應(yīng)力超過其屈服極限而破裂。破冰效果預(yù)測與所采用的破冰準(zhǔn)則有關(guān)，若采用文獻(xiàn)[12]提出的關(guān)系式作為破冰準(zhǔn)則，即：λ=max|?w/?x|≥0.04，則可以得到如圖6所示的破冰效果。此時(shí)可以看出，在氣墊船以亞、超臨界速度航行時(shí)破冰效果差，而以臨界速度航行時(shí)仍能得到與圖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 結(jié) 論

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

[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] 鹿飛飛, 張志宏, 胡明勇,等. 淺水岸壁條件下脈沖載荷引起的黏彈性浮冰層位移響應(yīng)[J]. 振動(dòng)與沖擊, 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] 張志宏, 鹿飛飛, 丁志勇,等. 勻速移動(dòng)載荷激勵(lì)浮冰層大幅響應(yīng)的臨界速度[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 盧再華, 張志宏, 胡明勇,等. 全墊升式氣墊船破冰過程的數(shù)值模擬[J]. 振動(dòng)與沖擊, 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] 劉巨斌, 張志宏, 張遼遠(yuǎn), 等. 氣墊船興波破冰問題的數(shù)值計(jì)算[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 劉巨斌, 張志宏, 張遼遠(yuǎn), 等. 邊界元—有限差分法在氣墊船破冰數(shù)值模擬中的應(yīng)用[J]. 海軍工程大學(xué)學(xué)報(bào), 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

國家自然科學(xué)基金(51479202)

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

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

張志宏 男,博士,教授,博士生導(dǎo)師,1964年生

O35

A

10.13465/j.cnki.jvs.2017.23.005