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

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

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

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

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

0　引　言

近年來，船舶朝著大型化、快速化迅速發(fā)展，隨著船舶噸位的增加，甲板面積的增大，船首外飄程度更大，由此帶來的首部砰擊問題更加嚴(yán)重。強(qiáng)烈的砰擊對(duì)船舶的危害很大，一方面使沖擊區(qū)域承受巨大的壓力，局部結(jié)構(gòu)屈曲甚至破壞；另一方面將引起整個(gè)船體劇烈的顫振，可能會(huì)導(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種不同的船首部的局部模型，研究其入水時(shí)的砰擊壓力峰值。文獻(xiàn)［5］建立了船首部剖面模型與船首部局部模型，對(duì)比了剖面垂直入水、船首部有航速與無航速入水 3 種情況的砰擊壓力峰值。文獻(xiàn)［6］建立了Aarsnes［7］以及 Arai&Matsunaga［8］進(jìn)行的船體剖面試驗(yàn)中的模型，并對(duì)其采用SPH方法進(jìn)行了入水砰擊計(jì)算。文獻(xiàn)［9］中對(duì)某大型集裝箱船首部局部模型的外飄砰擊壓力進(jìn)行了理論計(jì)算。以上文獻(xiàn)都沒有涉及全船計(jì)算模型的研究。

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

1　基本理論

1.1拉格朗日有限元法

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

改寫成

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

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

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

1.2歐拉有限體積法

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

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

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

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

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

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

2　仿真方法可靠性分析

為了驗(yàn)證 MSC.Dytran 軟件在計(jì)算流固耦合問題的可靠性，先對(duì) 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)格為水域，用無粘性、可壓縮的理想流體的材料填充。空氣域壓力用狀態(tài)方程 EOSGAM 描述，水區(qū)域壓力用多項(xiàng)式狀態(tài)方程 EOSPOL 描述。楔形體的外表面為流固耦合面，采用一般耦合算法，考慮到問題的二維性，模型在 z 方向上取單位長(zhǎng)度，且在歐拉單元的 z+，z- 前后表面上施加剛性墻邊界以保證流場(chǎng)僅在 xoy 平面內(nèi)運(yùn)動(dòng)，而楔形體也只允許有 y 方向上的運(yùn)動(dòng)。根據(jù) Chuang［1 1］的觀點(diǎn)，為了保證壓力波在流場(chǎng)中傳播不至于反射回來影響砰擊區(qū)域，本文所選取的水域的寬度滿足如下條件：

式中：WT為流場(chǎng)的寬度；Lm為楔形體的半寬。

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

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

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

圖 3　斜底最大砰擊壓力時(shí)歷曲線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ǔ)上，針對(duì)某大型集裝箱船進(jìn)行全船建模，計(jì)算船舶艏部舷側(cè)區(qū)域的外飄砰擊壓力。該船主尺度：總長(zhǎng)275.2 m，船寬 32.3 m，型深 21.8 m，設(shè)計(jì)吃水 12.2 m。

3.1計(jì)算工況

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

運(yùn)用基于非線性時(shí)域勢(shì)流理論的 Rankine 面元法，對(duì)集裝箱船型在上述工況下的運(yùn)動(dòng)響應(yīng)進(jìn)行預(yù)報(bào)。最終得到的船舶運(yùn)動(dòng)參數(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)動(dòng)參數(shù)Tab. 2　Ship motion parameters

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

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

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

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

4　結(jié)果分析

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

4.1航速對(duì)砰擊壓力極值大小的影響

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

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

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

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

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

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

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

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

4.2航速對(duì)砰擊壓力極值分布的影響

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

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

4.3原因分析

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

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

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

5　結(jié)　語

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

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

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

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

［1］戴仰山，沈進(jìn)威，宋競(jìng)正. 船舶波浪載荷［M］. 北京：國(guó)防工業(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)國(guó). 江海直達(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］劉正國(guó)，吳衛(wèi)國(guó)，潘晉，等. 基于船波相對(duì)運(yùn)動(dòng)的船首砰擊仿真方法［J］. 中國(guó)艦船研究，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］. 水動(dòng)力學(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動(dòng)力性能。