林 源，馬 寧，顧解忡

（1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240；2.高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海200240）

0 引言

隨著全球經(jīng)濟(jì)逐漸一體化，為實(shí)現(xiàn)航運(yùn)的規(guī)模效應(yīng)，集裝箱船的大型化日益加劇。集裝箱船尺度的增加，相應(yīng)地也會(huì)帶來(lái)甲板大開(kāi)口的現(xiàn)象，隨之帶來(lái)船體梁固有頻率下降，而航速的提升會(huì)增高波浪遭遇頻率，從而會(huì)引起船體的共振或強(qiáng)迫振動(dòng)，將影響船舶的使用壽命。所以在大型集裝箱船的設(shè)計(jì)制造過(guò)程中，水彈性響應(yīng)分析（如：波激振動(dòng)）應(yīng)當(dāng)被重視。

目前有關(guān)船舶水彈性理論已有較為全面的發(fā)展。最為經(jīng)典的水彈性理論就是Bishop 和Price[1]提出的二維理論，然而其忽略流體運(yùn)動(dòng)沿船長(zhǎng)方向的相互干擾，并未能計(jì)及船體首部和尾部的三維效應(yīng)[2]。接下來(lái)Wu[3]，Price 等[4]提出三維水彈性理論，是將三維耐波性與三維結(jié)構(gòu)動(dòng)力學(xué)相結(jié)合，但是三維理論需要大量的船體網(wǎng)格。

針對(duì)在工程上的快速水彈性預(yù)報(bào)，三維水彈性理論需要建立大量船體網(wǎng)格，工作量較大。因而簡(jiǎn)化網(wǎng)格對(duì)于數(shù)值計(jì)算是必然的，Golberg[5]和Fairweather 等[6]證明了MFS（The method of fundamental solutions）作為一種無(wú)網(wǎng)格方法在水動(dòng)力計(jì)算上的優(yōu)越性，但現(xiàn)階段MFS 方法尚未用于水彈性響應(yīng)的求解。而且有關(guān)的三維有限元與二維耐波性混合模型研究較少，Wang 等[7]將其應(yīng)用到浮體結(jié)構(gòu)中，但運(yùn)用到船舶響應(yīng)求解極少。本文提出一種混合水彈性模型，利用三維有限元方法預(yù)報(bào)船舶的固有模態(tài)，再基于MFS 方法的STF 切片理論[8-9]，進(jìn)行固有模態(tài)下的無(wú)網(wǎng)格船體廣義流體力計(jì)算，最后根據(jù)模態(tài)疊加法求解船舶在波浪作用下的波激振動(dòng)分析。為驗(yàn)證本方法的計(jì)算精度，采用Bernoulli 梁垂向振動(dòng)模態(tài)振型的解析解和水動(dòng)力軟件HydroSTAR 的數(shù)值計(jì)算結(jié)果進(jìn)行驗(yàn)證。

1 船體固有模態(tài)預(yù)報(bào)及波激振動(dòng)

根據(jù)船舶水彈性統(tǒng)一理論，得到船舶在波浪作用下振動(dòng)的響應(yīng)方程

其中：a，b 和c 分別表示船體廣義的結(jié)構(gòu)質(zhì)量、阻尼、剛度矩陣，是N×N 對(duì)角矩陣，每個(gè)矩陣分別由矩陣系數(shù)arr，brr和crr組成，通過(guò)結(jié)構(gòu)模態(tài)分析獲得。流體廣義質(zhì)量、阻尼、剛度矩陣分別由A、B 和C 表示，由矩陣系數(shù)Ars，Brs和Crs構(gòu)成，并通過(guò)基于MFS 方法的STF 切片方法進(jìn)行求解。pa和F 分別表示為主坐標(biāo)和廣義干擾力矩陣。ωe為遭遇頻率。

1.1 模態(tài)分析

1.1.1 干模態(tài)分析

干模態(tài)法是先把結(jié)構(gòu)系統(tǒng)與流體分開(kāi)，求出真空中結(jié)構(gòu)的固有頻率和固有模態(tài)，即干模態(tài)[10]。根據(jù)實(shí)際需要進(jìn)行模態(tài)分離，使自由度數(shù)目減少。基于干模態(tài)坐標(biāo)變換下，再進(jìn)行考慮流體的作用。

根據(jù)動(dòng)力平衡方程有

其中：M、C、K 分別為系統(tǒng)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣，F(xiàn)（t）為載荷向量，、、r 分別為節(jié)點(diǎn)的加速度、速度和位移向量。

真空中結(jié)構(gòu)無(wú)阻尼自由振動(dòng)方程如方程（3）所示，可通過(guò)求解得到相應(yīng)的固有頻率和固有模態(tài)。

假設(shè)系統(tǒng)做簡(jiǎn)諧運(yùn)動(dòng)，即為

其中：ω 為頻率，A 為振幅向量，將方程（4）代入方程（3）可得：

若使得方程（5）有非零解，A 不全為零，則有特征方程如下，解此方程得到的響應(yīng)的特征值和特征向量，即為相應(yīng)船體干模態(tài)的固有頻率ω 和固有振型。

MSC.Patran&Nastran 軟件[11]利用其經(jīng)驗(yàn)公式，已給出計(jì)算結(jié)果精度極高的方法。通過(guò)模態(tài)分析（Modal analysis），得到固有模態(tài)頻率ω 和固有模態(tài)振型，并導(dǎo)出振型歸一化的模態(tài)質(zhì)量陣a 和模態(tài)剛度陣c。對(duì)于結(jié)構(gòu)的廣義阻尼陣b 的相關(guān)計(jì)算可以按下列方程[10]進(jìn)行計(jì)算，其中

其中：v 為阻尼系數(shù)，ωr為r 階模態(tài)下的固有頻率，arr為r 階模態(tài)廣義質(zhì)量。小阻尼情況下，對(duì)數(shù)衰減率為δ≈2πv。一般情況下，其經(jīng)驗(yàn)公式如（8）式所示，從而可得到相應(yīng)的阻尼系數(shù)。

1.1.2 濕模態(tài)分析

船體的附連水振動(dòng)質(zhì)量在MSC.Patran&Nastran 內(nèi)是通過(guò)定義有限元模型濕表面單元號(hào)碼和船舶的吃水來(lái)自動(dòng)實(shí)現(xiàn)其計(jì)算，其理論是用Helmholtz 方法即源匯法（也叫邊界元法）解流體運(yùn)動(dòng)的拉普拉斯方程。濕模態(tài)是計(jì)及周?chē)畡?dòng)力影響的結(jié)構(gòu)系統(tǒng)的無(wú)阻尼自由振動(dòng)模態(tài)，可通過(guò)有限元軟件中設(shè)置相關(guān)的參數(shù)[11]求解濕模態(tài)。

1.2 廣義流體作用力

對(duì)于有航速的船體廣義流體質(zhì)量、阻尼和剛度矩陣元素[1]如下所示：

其中：m（x），N（x），B（x）分別為單位船體上的附加質(zhì)量、流體阻尼系數(shù)和微片寬度。對(duì)于廣義流體質(zhì)量、阻尼和剛度矩陣元素的積分，將Feng 等人[9]利用MFS 無(wú)網(wǎng)格方法求解船體水動(dòng)力系數(shù)的方法，應(yīng)用到基于固有模態(tài)下的廣義流體作用力系數(shù)的求解。

1.3 波激響應(yīng)

經(jīng)過(guò)上述計(jì)算，可以得到船體的廣義質(zhì)量陣a、阻尼陣b 和剛度陣c，及流體廣義質(zhì)量陣A、阻尼陣B、剛度矩陣C 和波浪激勵(lì)力陣F，便可列出船體的波激振動(dòng)模態(tài)分析方程[10]，即得到相應(yīng)的主坐標(biāo)Pa，從而可以求得船體在波浪作用下垂直彎曲振動(dòng)的動(dòng)位移w（x，t）、動(dòng)彎矩M（x，t）和動(dòng)切力V（x，t）等動(dòng)力響應(yīng)。

2 模型精度驗(yàn)證

為驗(yàn)證計(jì)算模型的精度，本文通過(guò)三維有限元固有模態(tài)計(jì)算結(jié)果與Bernoulli 梁的解析解[10]（如方程（15）所示）對(duì)比，從而對(duì)固有模態(tài)的求解精度進(jìn)行了驗(yàn)證。

圖1 Bernoulli 梁1 階模態(tài)的三維有限元計(jì)算結(jié)果與梁解析解對(duì)比Fig.1 The comparisons between the first order model shapes of Bernoulli beam by 3D FEM and the analytical solution

圖2 Bernoulli 梁2 階模態(tài)的三維有限元計(jì)算結(jié)果與梁解析解對(duì)比Fig.2 The comparisons between the second order model shapes of Bernoulli beam by 3D FEM and the analytical solution

圖3 Bernoulli 梁3 階模態(tài)的三維有限元計(jì)算結(jié)果與梁解析解對(duì)比Fig.3 The comparisons between the third order model shapes of Bernoulli beam by 3D FEM and the analytical solution

其中：kr為特征值，其取值滿足如下特征方程：cos krL·cosh krL=1。關(guān)于三維有限元結(jié)果和Bernoulli 梁的解析解，這里均采用了將船尾處的位移歸一化為單位1，方便對(duì)比分析。根據(jù)剖面上的點(diǎn)位移轉(zhuǎn)角和剖面存在一定的幾何關(guān)系，將三維Bernoulli 梁模型模態(tài)分析后剖面的節(jié)點(diǎn)的位移回歸得到歸一化的振型圖。圖1 至圖3 分別給出了1 至3 階固有模態(tài)計(jì)算結(jié)果和解析解的對(duì)比。其中1 階和2 階模態(tài)計(jì)算結(jié)果能夠與解析解達(dá)到90%以上的重合程度。從而證明應(yīng)用三維有限元方法進(jìn)行波激振動(dòng)的固有模態(tài)分析，較遷移矩陣隨迭代而逐步疊加誤差，能為進(jìn)一步基于模態(tài)振型下的流體水動(dòng)力系數(shù)求解提供更高的精度。

通過(guò)ANSYS 軟件對(duì)于三維有限元模態(tài)分析結(jié)果進(jìn)行交叉的計(jì)算驗(yàn)證，通過(guò)表1 可看出，MSC.Patran&Nastran計(jì)算結(jié)果有保證。

表1 MSC.Patran&Nastran 與ANSYS 計(jì)算固有模態(tài)對(duì)比分析Tab.1 Comparisons on modal analysis results between MSC.Patran &Nastran and ANSYS

圖4 大型集裝箱船垂蕩和縱搖附加質(zhì)量曲線Fig.4 Heave and pitch added mass of very large container ship

圖5 大型集裝箱船垂蕩和縱搖阻尼系數(shù)曲線Fig.5 Heave and pitch damping of very large container ship

本文為驗(yàn)證模型耐波性部分計(jì)算的正確性，計(jì)算了迎浪狀態(tài)下、100%服務(wù)航速時(shí)大型集裝箱船的垂蕩和縱搖的附加質(zhì)量、附加阻尼和激勵(lì)力，與商業(yè)軟件HydroSTAR 計(jì)算的結(jié)果對(duì)比分析，如圖4 至圖6所示。在中高頻波浪下，本文的計(jì)算結(jié)果精度可以滿足，對(duì)于由于切片法本身的弊端，低頻情況下的精度上仍需改善。但是總體上的誤差在15%以?xún)?nèi)，并且由于低頻部分對(duì)于波激振動(dòng)趨勢(shì)的展示影響較小，故而可以滿足工程需求。

圖6 大型集裝箱船垂蕩和縱搖波浪激勵(lì)力曲線Fig.6 Heave and pitch excitation of very large container ship

3 大型集裝箱船波激振動(dòng)分析

對(duì)大型集裝箱船的有限元分析，是通過(guò)MSC.Patran &Nastran 軟件，得到了各階固有干模態(tài)的固有振型及相應(yīng)的廣義質(zhì)量、廣義剛度，并通過(guò)經(jīng)驗(yàn)公式得到相應(yīng)的廣義阻尼，匯總成表，如表2所示。

表2 大型集裝箱船的干濕模態(tài)預(yù)報(bào)結(jié)果及結(jié)構(gòu)廣義質(zhì)量、剛度、阻尼值Tab.2 Dry and wet mode prediction and the generalized mass,damping and stiffness

根據(jù)剖面上的點(diǎn)位移轉(zhuǎn)角和剖面存在一定的幾何關(guān)系，將三維集裝箱船模型模態(tài)分析后剖面的節(jié)點(diǎn)的位移回歸得到歸一化的振型圖。其中大型集裝箱船的第1 階至第3 階對(duì)稱(chēng)振動(dòng)的固有振型圖，如圖7 至圖9所示，考慮到大型集裝箱船的寬大型、大甲板開(kāi)口復(fù)雜性，采用三維有限元方法對(duì)船舶進(jìn)行模態(tài)分析，計(jì)算結(jié)果與船體梁的遷移矩陣模態(tài)分析相比，精準(zhǔn)度較高。

圖7 某大型集裝箱船1 階垂向模態(tài)振型圖預(yù)報(bào)Fig.7 The predicted results on first order mode shape of very large container ship

圖8 某大型集裝箱船2 階垂向模態(tài)振型圖預(yù)報(bào)Fig.8 The predicted results on second order mode shape of very large container ship

圖9 某大型集裝箱船3 階垂向模態(tài)振型圖預(yù)報(bào)Fig.9 The predicted results on third order mode shape of very large container ship

通過(guò)上述模態(tài)分析，將回歸后的歸一化振型運(yùn)用到廣義流體作用力的計(jì)算，從而進(jìn)一步可以求解相應(yīng)的波激對(duì)稱(chēng)振動(dòng)方程。為考慮航速對(duì)于主坐標(biāo)的影響，采用迎浪下100%服務(wù)航速和50%服務(wù)航速的數(shù)值計(jì)算，可以通過(guò)圖10 看到對(duì)稱(chēng)響應(yīng)主坐標(biāo)幅值隨船長(zhǎng)與波長(zhǎng)之比L/λ 變化趨勢(shì)。對(duì)于第1階和第2 階垂直彎曲模態(tài)，其主坐標(biāo)變化趨勢(shì)相同，且第1 階主坐標(biāo)最大幅值大于第2 階主坐標(biāo)最大幅值，具有船舶波激響應(yīng)的一般趨勢(shì)。當(dāng)航速增大時(shí)，主坐標(biāo)的數(shù)值相應(yīng)也有一定程度的增大，并且第一個(gè)峰值皆出現(xiàn)在L/λ=1 處，證明了計(jì)算結(jié)果的準(zhǔn)確性。隨著遭遇頻率的增加，在L/λ≈8.75 和18.68 時(shí)產(chǎn)生第二個(gè)和第三個(gè)峰值。根據(jù)波長(zhǎng)與波頻關(guān)系式有L/λ=Lω2/2πg(shù)，遭遇頻率ωe=ω-ω2Ucos β/g，可以求出對(duì)應(yīng)L/λ≈8.75、18.68 時(shí)，遭遇頻率為2.353 rad/s 和4.147 rad/s，說(shuō)明本方法的可行性。而且較三維水彈性方法，采用無(wú)網(wǎng)格MFS 的STF 法求解流體水動(dòng)力系數(shù)可簡(jiǎn)化網(wǎng)格，降低計(jì)算難度和節(jié)省時(shí)間。

圖10 不同航速迎浪航行時(shí)的主坐標(biāo)和變化曲線（左圖為，右圖為）Fig.10 The third and fourth principal mode curves of various speeds in the heading sea

4 結(jié)論

本文以集裝箱船為目標(biāo)船型，提出了三維有限元結(jié)構(gòu)分析與二維耐波性分析相結(jié)合的水彈性預(yù)報(bào)方法，并基于該方法對(duì)某大型集裝箱船的波激響應(yīng)進(jìn)行了分析：

（1）本文方法較目前的三維水彈性分析大幅簡(jiǎn)化了網(wǎng)格劃分、提高了計(jì)算效率；且與二維水彈性方法相比，在結(jié)構(gòu)分析時(shí)能夠提高精度，具有較好的工程實(shí)用價(jià)值；

（2）本文通過(guò)MSC.Patran &Nastran 軟件模態(tài)分析得到的固有振型方法，經(jīng)過(guò)梁模型驗(yàn)證，平均精度在90%以上；計(jì)算的固有頻率，與ANSYS 軟件計(jì)算結(jié)果相比，精度有保證；

（3）本文基于STF 切片理論計(jì)算的廣義流體相關(guān)系數(shù)，忽略船體變形的剛體運(yùn)動(dòng)，與三維HydroSTAR軟件對(duì)比，精度在85%以上，具有工程實(shí)用價(jià)值，但在未來(lái)階段仍需改進(jìn)低頻部分的精度問(wèn)題；

（4）基于本文模型得到的波激對(duì)稱(chēng)響應(yīng)，符合船舶水彈性相關(guān)特性，得到共振的遭遇頻率分別為2.353 rad/s 和4.147 rad/s。