肖考考，高霄鵬（海軍工程大學(xué) 艦船工程系，湖北　武漢 430033）

基于全船建模的航速對船首外飄砰擊影響研究

肖考考，高霄鵬
（海軍工程大學(xué) 艦船工程系，湖北武漢 430033）

建立二維楔形體入水模型，驗(yàn)證入水砰擊仿真方法可靠性，對集裝箱船體進(jìn)行全船體建模，導(dǎo)入船體運(yùn)動參數(shù)，增加船體的縱搖運(yùn)動，采用一般耦合算法（General coupling），流體域用 Euler 單元模擬，船體設(shè)為剛體，劃分為 Lagrange 有限元網(wǎng)格，對船體進(jìn)行入水仿真。對比不同工況下的結(jié)果表明：全船模擬時船型縱向斜升角會發(fā)生變化，導(dǎo)致航速對船首入水砰擊壓力的影響較大，隨著航速的增加，入水砰擊壓力變大，同時航速還會使砰擊壓力峰值位置發(fā)生變化。

外飄砰擊；全船建模；一般耦合法；縱向斜升角

0　引　言

近年來，船舶朝著大型化、快速化迅速發(fā)展，隨著船舶噸位的增加，甲板面積的增大，船首外飄程度更大，由此帶來的首部砰擊問題更加嚴(yán)重。強(qiáng)烈的砰擊對船舶的危害很大，一方面使沖擊區(qū)域承受巨大的壓力，局部結(jié)構(gòu)屈曲甚至破壞；另一方面將引起整個船體劇烈的顫振，可能會導(dǎo)致船舶總縱強(qiáng)度或者局部強(qiáng)度的喪失［1］。

目前，關(guān)于船舶首部砰擊的研究，多采用首部局部建?；蛘叽灼拭娼?，文獻(xiàn)［2］采用的建立首部局部模型研究 S175 型船首部甲板上浪問題。文獻(xiàn)［3］運(yùn)用MSC.Dytran 模擬船體二維剖面入水來計(jì)算船體底部砰擊壓力，沒有涉及三維船體的計(jì)算。文獻(xiàn)［4］建立了 3種不同的船首部的局部模型，研究其入水時的砰擊壓力峰值。文獻(xiàn)［5］建立了船首部剖面模型與船首部局部模型，對比了剖面垂直入水、船首部有航速與無航速入水 3 種情況的砰擊壓力峰值。文獻(xiàn)［6］建立了Aarsnes［7］以及 Arai&Matsunaga［8］進(jìn)行的船體剖面試驗(yàn)中的模型，并對其采用SPH方法進(jìn)行了入水砰擊計(jì)算。文獻(xiàn)［9］中對某大型集裝箱船首部局部模型的外飄砰擊壓力進(jìn)行了理論計(jì)算。以上文獻(xiàn)都沒有涉及全船計(jì)算模型的研究。

全船建模相比于船體剖面建模和船首部局部建模的優(yōu)勢在于能更準(zhǔn)確地預(yù)報其船體運(yùn)動響應(yīng)，并可將預(yù)報結(jié)果直接作為砰擊預(yù)報中的船體運(yùn)動參數(shù)輸入；其次，在仿真中可以附加船舶的縱搖運(yùn)動，綜合體現(xiàn)船體在波浪中縱搖與垂蕩運(yùn)動的耦合運(yùn)動狀態(tài)，更接近實(shí)船在波浪中的運(yùn)動模式。故本文采用全船建模的方式對航速對船舶入水砰擊的影響問題進(jìn)行研究。

1　基本理論

1.1拉格朗日有限元法

MSC.Dytran 是瞬態(tài)動力學(xué)流固耦合領(lǐng)域的高端軟件，它的求解方法在時間域上采用顯式時間積分法。設(shè)當(dāng)前時間步是 n，顯式積分方法是將運(yùn)動微分方程

改寫成

如果將單元質(zhì)量分布在節(jié)點(diǎn)上，則為一對角陣，稱為集中質(zhì)量矩陣，線性方程組將成為一系列關(guān)于各個自由度的獨(dú)立的一元一次方程，從而可求出節(jié)點(diǎn)加速度：

如果假設(shè)加速度在一個時間步長內(nèi)恒定，在時間推進(jìn)上采用中心差分法：

顯式積分法不需要做矩陣分解，因此具有很高的計(jì)算效率。

1.2歐拉有限體積法

初始條件已知情況下，利用控制方程求解每一歐拉單元在一個時間步的密度、速度、比內(nèi)能和壓力。假設(shè) tn時刻各物理參數(shù)已知，對相鄰單元形心處流速進(jìn)行線性插值求出單元邊界處流速

求出穿越單元表面的質(zhì)量、動量及能量的流量

式中：ρ2為相鄰單元密度；et為單元質(zhì)量的總能量；T為邊界上單元面積上的面力；△V 為時刻 tn～tn + 1的一個時間步長內(nèi)穿越該單元的表面體積流量。采用單點(diǎn)高斯積分，通過控制方程從而可以解出單元形心處的物理量（密度、流速、內(nèi)能）在 tn + 1時刻的值，得到從時刻 tn至 tn + 1的變化量關(guān)系。根據(jù)材料本構(gòu)關(guān)系，可以進(jìn)一步計(jì)算出壓力值。

1.3拉格朗日-歐拉流固耦合

MSC.Dytran 中的流固耦合計(jì)算就是拉格朗日域（固體）與歐拉域（流場）的耦合計(jì)算。本文采用的一般耦合法大多是拉格朗日的固體在歐拉的流場范圍內(nèi)運(yùn)動，即拉格朗日域驅(qū)動歐拉域。流場雖有速度，但在流固耦合過程中歐拉網(wǎng)格不移動也不變形。

耦合的目的是為了讓歐拉網(wǎng)格與拉格朗日網(wǎng)格發(fā)生相互作用。在默認(rèn)設(shè)定下，歐拉求解器與拉格朗日求解器單獨(dú)分開。當(dāng)沒有定義流體和固體之間耦合關(guān)系時，即使拉格朗日單元恰好處在歐拉網(wǎng)格范圍內(nèi)也不會對歐拉材料的流動產(chǎn)生任何影響，同時自身也不受任何來自歐拉材料力的作用。耦合關(guān)系的定義可以啟動歐拉與拉格朗日單元之間的相互作用的計(jì)算程序，從而能夠分析復(fù)雜的流固禍合問題。

2　仿真方法可靠性分析

為了驗(yàn)證 MSC.Dytran 軟件在計(jì)算流固耦合問題的可靠性，先對 N.Aquelet［10］二維剛性楔形體（斜升角B = 10°）等速入水模型進(jìn)行仿真計(jì)算。本文建立了完整的二維模型，如圖 1 所示。

圖 1　計(jì)算模型Fig. 1　Calculation model

楔形體結(jié)構(gòu)用 Lagrange 四邊形殼單元來離散，采用剛性材料，以 5.425 m /s 的速度等速入水，底部距液面的初始高度為 0.02 m，流體材料用六面體的 Euler 單元來離散。為提高計(jì)算效率，在離楔形體較近的位置歐拉網(wǎng)格劃分較密，離楔形體較遠(yuǎn)的位置網(wǎng)格劃分較為稀疏。Euler，Lagrange 單元的尺寸比約為 1：1.2。如圖 1 所示，上層歐拉網(wǎng)格為空氣域，用可壓縮理想氣體的材料來填充；下層歐拉網(wǎng)格為水域，用無粘性、可壓縮的理想流體的材料填充?？諝庥驂毫τ脿顟B(tài)方程 EOSGAM 描述，水區(qū)域壓力用多項(xiàng)式狀態(tài)方程 EOSPOL 描述。楔形體的外表面為流固耦合面，采用一般耦合算法，考慮到問題的二維性，模型在 z 方向上取單位長度，且在歐拉單元的 z+，z- 前后表面上施加剛性墻邊界以保證流場僅在 xoy 平面內(nèi)運(yùn)動，而楔形體也只允許有 y 方向上的運(yùn)動。根據(jù) Chuang［1 1］的觀點(diǎn)，為了保證壓力波在流場中傳播不至于反射回來影響砰擊區(qū)域，本文所選取的水域的寬度滿足如下條件：

式中：WT為流場的寬度；Lm為楔形體的半寬。

圖 2 為楔形體入水時水域壓力云圖；圖 3 為斜底最大砰擊壓力的時間歷程，圖中曲線（1）為本文仿真解，曲線（2）為文獻(xiàn)［10］理論解，曲線（3）為文獻(xiàn)［10］數(shù)值解；圖 4 為各個階段自由液面的變化情況。

讀取圖中時間歷程數(shù)據(jù)可得，本文仿真解最大砰擊壓力為 1.286 MPa，而文獻(xiàn)中的理論解的最大砰擊壓力為 1.2 MPa，數(shù)值解最大砰擊壓力為 1.22 MPa，本文仿真解略大于文獻(xiàn)中的解，相對于理論解的誤差為7.16%，相對于數(shù)值解的誤差為 5.40%，均在 7%左右，說明本文仿真計(jì)算具有可靠性。

圖 2　水域的壓力云圖Fig. 2　Pressure cloud image

圖 3　斜底最大砰擊壓力時歷曲線Fig. 3　The curve of the maximum pressure

圖 4　楔形體入水自由液面變化情況Fig. 4　The change of the water level of the wedge into the water

3　船舶砰擊計(jì)算模型

在 MSC.Dytran 軟件計(jì)算入水砰擊問題可靠性的基礎(chǔ)上，針對某大型集裝箱船進(jìn)行全船建模，計(jì)算船舶艏部舷側(cè)區(qū)域的外飄砰擊壓力。該船主尺度：總長275.2 m，船寬 32.3 m，型深 21.8 m，設(shè)計(jì)吃水 12.2 m。

3.1計(jì)算工況

針對船舶航行過程中較危險海況及常見航速，選取了 6 種工況進(jìn)行計(jì)算，計(jì)算工況如表 1 所示。

運(yùn)用基于非線性時域勢流理論的 Rankine 面元法，對集裝箱船型在上述工況下的運(yùn)動響應(yīng)進(jìn)行預(yù)報。最終得到的船舶運(yùn)動參數(shù)見表 2。

3.2參數(shù)設(shè)置

建立了整船模型，如圖 5 所示，并將船體模型進(jìn)行網(wǎng)格劃分，建立水域并進(jìn)行計(jì)算參數(shù)的設(shè)置。

表 1　計(jì)算工況Tab. 1　Calculation condition

表 2　船舶運(yùn)動參數(shù)Tab. 2　Ship motion parameters

圖 5　全船模型Fig. 5　Whole ship model

船體結(jié)構(gòu)采用拉格朗日 4 節(jié)點(diǎn)四邊形單元劃分。船體網(wǎng)格劃分模型如圖 6 所示。由于不考慮結(jié)構(gòu)的局部變形和結(jié)構(gòu)在砰擊作用下在水中的振動，故將船體結(jié)構(gòu)設(shè)置為剛體材料。水域單元為歐拉六面體實(shí)體單元。如圖 6 所示，上層歐拉網(wǎng)格為空氣域，下層歐拉網(wǎng)格為水域。Euler 與 Lagrange 單元的尺寸比為 0.6：1。

模型入水運(yùn)動狀態(tài)定義為強(qiáng)迫運(yùn)動，模型運(yùn)動由船舶沿航向的運(yùn)動、垂蕩運(yùn)動以及縱搖運(yùn)動耦合而成，更真實(shí)的模擬了船舶在實(shí)際航行過程中發(fā)生砰擊時的運(yùn)動狀態(tài)，同時也考慮了全船仿真時入水砰擊的三維特性。將表 2 中的運(yùn)動參數(shù)作為船舶運(yùn)動的輸入條件。

圖 6　有限元模型Fig. 6　Finite element model

4　結(jié)果分析

如圖 7 所示，取船舷側(cè) 0～3 站，16～19 號水線上的16 個預(yù)報點(diǎn)讀取砰擊壓力值。下面列出有義波高 14 m，航速分別為 0 kn，10 kn，18 kn 三種工況下某預(yù)報點(diǎn)的砰擊壓力時歷曲線，如圖 8～圖 10 所示。

4.1航速對砰擊壓力極值大小的影響

在圖 8～圖 10 中，對比有義波高為 14 m 時，航速18 kn，10 kn，0 kn 兩種工況的砰擊壓力值，可以看出，同樣海況下，有航速時預(yù)報點(diǎn)的砰擊壓力值普遍大于無航速情況。

圖 7　外飄砰擊預(yù)報點(diǎn)位置示意圖Fig. 7　A sketch map of the prediction point of the flare slamming

圖 8　航速 0 kn 下 2-1 點(diǎn)時歷曲線Fig. 8　2-1 point time history curve at speed of 0kn

圖 9　航速 10 kn 時 2-1 點(diǎn)時歷曲線Fig. 9　2-1 point time history curve at speed of 10kn

圖 10　航速 18 kn 時 2-1 點(diǎn)時歷曲線Fig. 10　2-1 point time history curve at speed of 18kn

同樣，由圖 11 可看出，航速 18 kn，有義波高 14 m工況下的砰擊壓力極值比航速 10 kn，有義波高 18.5 m工況下砰擊壓力極值要高。一般情況下，海況越高，船舶受到的砰擊壓力越大，而此時當(dāng)海況由浪高 18.5 m降至 14 m 時，航速提高至 18 kn 反而使砰擊壓力更大，說明航速是影響船舶砰擊壓力大小的重要因素之一。

圖 11　極值沿舷側(cè)高度分布曲線Fig. 11　Extreme value side along the height distribution curve

當(dāng)考慮船體的縱搖運(yùn)動時，航速對于船首砰擊壓力的影響顯著，可總結(jié)其規(guī)律為：航速越高，船首砰擊壓力越大。

4.2航速對砰擊壓力極值分布的影響

從圖 12 可以觀察到，總共 6 個不同工況下的砰擊壓力極值曲線，6 個不同工況包含了 4 個不同航速，可以發(fā)現(xiàn)，航速為 0 kn 與 10 kn 的曲線的峰值發(fā)生在3-1 預(yù)報點(diǎn)處，航速為 18 kn 的曲線的峰值發(fā)生在 2-1預(yù)報點(diǎn)處，航速為 26 kn 的曲線的峰值發(fā)生在 1-1 預(yù)報點(diǎn)處?？梢园l(fā)現(xiàn)，航速越高，砰擊壓力極值發(fā)生的點(diǎn)越往船尾推移。且航速為 0 kn 的曲線，砰擊壓力值增長的速度比航速為 10 kn 的曲線要慢，可以推斷，航速 0 kn 下比航速 10 kn 下砰擊壓力峰值發(fā)生點(diǎn)更靠近船首。由此可以總結(jié)得出航速對于船首砰擊壓力極值分布的影響為：航速越高，船首砰擊壓力極值發(fā)生點(diǎn)越向船尾方向推移。由于圖 12 曲線中，航速相同，海況不同時峰值發(fā)生位置相同，故可進(jìn)一步推斷，砰擊壓力極值點(diǎn)發(fā)生位置在船長方向上并不受海況的影響，只與航速的高低有關(guān)。

圖 12　極值沿船長方向分布曲線Fig. 12　Extreme value side along the direction of ship distribution curve

4.3原因分析

分析其主要原因是，全船建模仿真時，將船舶的縱搖運(yùn)動考慮進(jìn)去之后，入水運(yùn)動會改變船舶首部區(qū)域縱向剖面與水面的縱向斜升角，當(dāng)船首入水時，會隨著船舶縱搖角度的增大而減小，當(dāng)縱搖角速度一定時，航速越大，船舶縱向剖面的入水速度也越大，由于縱向斜升角的減小，根據(jù) Wagner［12］的理論，入水沖擊壓力寫為：

式中：ρ 為水密度；V 為入水速度；K 為無因次壓力系數(shù)，取值根據(jù)有效沖擊角的范圍而定。

由文獻(xiàn)［11］中總結(jié)的結(jié)論，物體入水角度變小，則砰擊壓力系數(shù) K 增大。當(dāng) K 與 V 同時增大時，入水沖擊壓力 pi會顯著增大，則縱向剖面入水產(chǎn)生的壓力會明顯變大。

5　結(jié)　語

本文基于船舶全船建模，對某集裝箱船外飄砰擊過程進(jìn)行了數(shù)值仿真，得到如下結(jié)論：

1）全船入水仿真時，縱搖運(yùn)動導(dǎo)致船型縱向斜升角變小，從而導(dǎo)致航速對于船首外飄砰擊壓力大小的影響較大，其規(guī)律為航速越高，船首砰擊壓力越大。

2）全船入水仿真時，航速越高，船首砰擊壓力極值發(fā)生點(diǎn)越向船尾推移。

3）與船首局部建?；蚨S船體剖面建模相比，本文中嘗試的將船舶的航速與縱搖運(yùn)動考慮進(jìn)去的全船建模仿真方法能更貼近實(shí)船運(yùn)動模式，考慮的因素更多，更全面。

［1］戴仰山，沈進(jìn)威，宋競正. 船舶波浪載荷［M］. 北京：國防工業(yè)出版社，2007.

［2］向紅貴. 高速船甲板上浪的數(shù)值模擬研究［D］. 上海：上海交通大學(xué)，2008. XIANG Hong-gui. Numerical simulation of green water occurrence of high speed vessels［D］. Shanghai：Shanghai Jiaotong University，2008.

［3］文志飛. 船舶砰擊載荷的計(jì)算方法研究［D］. 哈爾濱：哈爾濱工程大學(xué)，2010. WEN Zhi-fei. A computational method for slamming loads of ship［D］. Harbin：Harbin Engineering University，2010.

［4］孫峰，吳衛(wèi)國. 江海直達(dá)船船首入水砰擊研究［J］. 交通科技，2011（5）：114-116. SUN Feng，WU Wei-guo. Impacting study on the water entry of the bow of river-sea ship［J］. Transportation Science & Technology，2011（5）：114-116.

［5］劉正國，吳衛(wèi)國，潘晉，等. 基于船波相對運(yùn)動的船首砰擊仿真方法［J］. 中國艦船研究，2013，8（6）：20-26. LIU Zheng-guo，WU Wei-guo，PAN Jin，et al. Simulation method for the slamming problems of ship bows based on the relative motion between hulls and wave surfaces［J］. Chinese Journal of Ship Research，2013，8（6）：20-26.

［6］胥飛，劉樺. 艏外飄型船體入水砰擊的二維SPH模擬［J］. 水動力學(xué)研究與進(jìn)展，2013，28（5）：585-590. XU Fei，LIU Hua. Numerical simulation of 2-D bow flare slamming using SPH method［J］. Chinese Journal of Hydrodynamics，2013，8（5）：585-590.

［7］AARSNES J V. Drop test with ship sections-effect of roll angle［R］. Marintek Report，Trondheim：Norwegian Marine Technology Research Institute，1996.

［8］ARAI M，TASAKI R. A numerical study of water entrance of two-dimensional wedges-effect of gravity，spray generation and vertical load［C］//Proceedings of the 3rd International Symposium on Practical Design of Ships and Mobile Units. Trondheim，Norway：Norwegian Institute of Technology，1987.

［9］陳震，馮永軍，肖熙. 大型集裝箱船舷側(cè)外飄砰擊特性研究［J］. 船海工程，2011，40（3）：1-4，9. CHEN Zhen，F(xiàn)ENG Yong-jun，XIAO Xi. Study on flare impact characteristics of large container ships［J］. Ship & Ocean Engineering，2011，40（3）：1-4.

［10］AQUELET N，SOULI M，OLOVSSON L. Euler-Lagrange coupling with damping effects：application to slamming problems［J］. Computer Methods in Applied Mechanics and Engineering，2006，195（1/3）：110-132.

［11］STAVOVY A B，CHUANG S L. Analytical determination of slamming pressures for high-speed vehicles in waves［J］. Journal of Ship Research，1967，20（4）：190-198.

［12］Wagner H. über Sto?- und Gleitvorg?nge an der Oberfl?che von Flüssigkeiten［J］. Zeitschrift für Angewandte Mathematik und Mechanik，1932，12（4）：193-215.

Effect of ship speed for the flare slamming of ship bows based on the whole ship modeling

XIAO Kao-kao，GAO Xiao-peng
（Naval University of Engineering，Department. of Naval Architecture & Ocean Engineering，Wuhan 430033，China）

The impact of one rigid wedge subject is calculated to verify the reliability of the slamming simulation method.Modeling a whole ship，import the ship's motion parameters，including pitching motion. In simulation，the fluid is represented by an Euler formulation with 8-nodes brick elements，the ship structure is represented by a Lagrange grid，and simulate the whole ship model's slamming problems by using General coupling. The results show that：In case of high speed，moderate sea conditions，the bow's slamming pressure is higher than the other conditions；the longitudinal deadrise angle of the bow will be changed in whole model simulation，resulting in a greater influences the ship speed to the bow's slamming pressures，and the higher of ship speed，the higher of the bow's slamming pressure and the ship speed also makes slamming pressure peak position changes.

flare slamming；whole ship modeling；general coupling；longitudinal deadrise angle

U661

A

1672 - 7619（2016）08 - 0023 - 06

10.3404/j.issn.1672 - 7619.2016.08.005

2015 - 09 - 17；

2015 - 10 - 08

航空基金資助項(xiàng)目（2013102221284811）

肖考考（1991 - ），男，碩士研究生，研究方向?yàn)榕灤黧w動力性能。