翟鋼軍， 程 勇， 馬 哲

(大連理工大學(xué) 深海工程研究中心,遼寧 大連 116024)

隨著我國工業(yè)化進程的加快和國民經(jīng)濟的迅速發(fā)展，我國的石油消費量呈迅速上升的勢頭，而我國石油資源供應(yīng)能力有限，因此借鑒國際先進經(jīng)驗，迅速建立國家石油戰(zhàn)略儲備，是最現(xiàn)實、最可靠的手段。目前石油儲備主要是針對原油儲備，其方式有5種[1]：地上浮頂油罐儲備方式、地中油罐儲備方式、地下巖洞儲備方式、地下鹽穴儲備方式和海上石油儲備方式。結(jié)合我國實際情況，采用以海上戰(zhàn)略石油儲備為主的多樣化儲備方式更加適合我國人口密度、地理環(huán)境和工業(yè)分布狀況的特點。

圖1 海上石油儲備基地布置圖

海上石油儲備基地如圖1所示，主要由數(shù)艘儲油船并排靠泊組成，其中系泊系統(tǒng)由若干系纜樁和橡膠護舷構(gòu)成。儲油船屬于超大型浮式結(jié)構(gòu)物(VLFS)，即水平尺度與高度的比值非常大，是一個極為扁平的柔性結(jié)構(gòu)物[2]。由于儲油船的彈性變形與剛體位移幾乎屬于同一量級，計算結(jié)構(gòu)在流體作用下的水動力特性時，不能再把船體當(dāng)作單一的剛體進行計算，必須分析儲油船在流體荷載作用下的底部變形，以及船體變形對流場的影響，這是一個典型的流固耦合問題，需用水彈性理論進一步分析。

水彈性力學(xué)成立于70年代，由Bishop等[3]提出，它需要考慮慣性力以及流體荷載和結(jié)構(gòu)彈性力之間的耦合作用，關(guān)系到水動力學(xué)和結(jié)構(gòu)動力學(xué)兩門學(xué)科以及它們之間的聯(lián)系，屬于典型的交叉學(xué)科。日本、美國、英國和韓國等相繼對超大型浮式結(jié)構(gòu)物的水彈性進行了研究，推行了一系列VLFS研究計劃，如Wu等[4]通過特征函數(shù)展開法來求解彈性薄板的二維水彈性問題；Hong等[5]分析了連接OWC防波堤的VLFS的頻域水彈性響應(yīng)和漂移力；Murai等[6]提出了一種子結(jié)構(gòu)法來加速水彈性計算且不丟失原有精度；Fu等[7]預(yù)測了兩個柔性連接的VLFS水彈性響應(yīng)；Watanabe 等[8]提出了一個圓柱形VLFS的標準水彈性方法；Pham 等[9]進一步分析了配有減振裝置的圓柱形VLFS的水彈性問題。我國也有很多學(xué)者從事該方面研究，如吳有生[10]發(fā)展了適用非常廣泛的三維水彈性理論；崔維成等[2,11-12]總結(jié)了VLFS的研究現(xiàn)狀和一些關(guān)鍵問題，并運用模態(tài)展開法分別計算了梁模型和矩形平板模型的水彈性響應(yīng)；滕斌等[13]學(xué)者介紹了計算水彈性的模態(tài)展開法和直接計算法兩種方法，并計算了五種振型函數(shù)下結(jié)構(gòu)的水彈性動力學(xué)；張淑華等[14]學(xué)者采用直接法計算了彈性平板模型的撓度幅值。

對于超大型浮式儲油船的水動力分析，上海交通大學(xué)的李文龍等[16]基于三維勢流理論，建立了儲油船的剛體運動方程，并計算了相關(guān)的水動力參數(shù)。然而儲油船屬于超大型結(jié)構(gòu)，其彈性變形甚至比剛性位移更為重要，因此必須給予考慮。本文采用模態(tài)分析法，將船的水彈性位移表示成一系列模態(tài)振型的疊加，應(yīng)用邊界元法求出繞射勢和一系列模態(tài)上單位速度運動引起的輻射勢。再通過儲油船的彈性體運動方程求出各模態(tài)的振型幅值，最后回代到水彈性位移表達式中，求出儲油船的水彈性響應(yīng)。在儲油船的水彈性力學(xué)分析中，本文假設(shè)儲油船為一薄板模型，通過與滕斌等[13]的數(shù)值結(jié)果及Maeda等[15]的實驗結(jié)果進行對比，驗證了本文方法的可行性，進一步分析了不同水深及波要素對船體運動響應(yīng)及水動力參數(shù)的影響。

1 儲油船的水彈性分析

考慮波浪與儲油船相互作用時，儲油船長為L，寬為B，吃水為d，且開敞水域水深為h，選定坐標系如圖2所示，x-y平面位于未干擾自由水面，根據(jù)右手定則z軸垂直向上。整個流域Ω由儲油船底面Sb，側(cè)面Sc，自由水面Sf，無窮遠柱面S∞和海底面Sd面組成。

儲油船在流場中受流體作用過程中，包括結(jié)構(gòu)和流場兩部分計算，船受到流體荷載作用要產(chǎn)生變形，屬于結(jié)構(gòu)計算，而船體的存在和變形改變了流場，屬于流場計算，反之，流場的變化也對船體的彈性變形產(chǎn)生影響，總之兩者相互聯(lián)系，相互耦合。其中對流域的求解可采用邊界元法，對結(jié)構(gòu)的彈性運動方程求解可采用直接法或模態(tài)分析法(本文采用該方法)，下面將對這兩個方面進行詳細論述。

圖2 坐標示意圖

1.1 流場控制方程和邊界條件

假定流體為無旋、不可壓縮的理想流體，則流體運動可以通過速度勢Φ來表述。當(dāng)考慮入射波浪為一頻率為ω的規(guī)則波時，所有含時間變量t的物理量均可分離出時間分量，表達如下：

Φ=Re{φ(x,y,z)e-iωt}

(1a)

W(x,y,t)=Re{w(x,y)e-iωt}

(1b)

P(x,y,t)=Re{p(x,y)e-iωt}

(1c)

式中：W(x,y,t)為船的垂向運動響應(yīng)；P(x,y,t)為船底部所受的流體壓力。

進一步假定波浪的波陡和船體運動量足夠小，滿足線性勢流理論。因此，復(fù)速度勢φ滿足以下邊界值定解問題：

2φ=0 在整個流域

(2)

(3)

(4)

(5)

(6)

ω2=gktanh(kh)

(7)

式中：n為單位法向量，指出流域為正；r為參考點到坐標原點的水平距離；復(fù)值變量w(x,y)為儲油船底部的垂直位移，包括剛體位移和彈性變形兩部分。

1.2 儲油船運動方程

對于超大型浮式儲油船，水平尺度遠大于垂直尺度，因此可將其簡化為彈性薄板，則問題轉(zhuǎn)化為求解彈性薄板在波浪作用下的水彈性響應(yīng)。典型的薄板振動微分方程為：

D4W(x,y,t)=Q(x,y,t)

(8)

式中：D為薄板的彎曲剛度；Q為薄板單位面積上的橫向荷載，包括慣性力和流體作用力兩部分，表達如下：

(9)

式中：ms為薄板單位面積質(zhì)量，W(x,y,t)和P(x,y,t)表達如上。將式(1b),式(1c)代入到式(8)中可得儲油船的彈性運動方程為：

D4w(x,y)-ω2msw(x,y)=p(x,y)

(10)

對于自由漂浮的儲油船，其端部邊界上的彎矩和剪力為零：

(11)

(12)

式中：v是泊松比。

求出彈性位移后，儲油船內(nèi)部彎矩可求得為：

(13a)

(13b)

2 儲油船的運動響應(yīng)及速度勢的模態(tài)展開形式

由儲油船的彈性運動方程(10)可以看出，結(jié)構(gòu)運動響應(yīng)w(x,y)和流體運動是耦合在一起的，并且由物面條件式(5)可知流體運動可由結(jié)構(gòu)運動響應(yīng)得到。因此，將問題分解為速度勢的水動力學(xué)問題和薄板振動的結(jié)構(gòu)動力學(xué)問題，船的運動響應(yīng)和速度勢都進行模態(tài)展開。

2.1 振型函數(shù)的一般形式

對于儲油船的振型函數(shù)可以表示成x和y兩方向振型函數(shù)的乘積形式，因此以一維梁模型為例。振型函數(shù)的選取分為干模態(tài)法和濕模態(tài)法，本文采取前者方法。干模態(tài)法中，以在空氣中自由振動的振動模態(tài)為自然模態(tài)。在水彈性問題中，自然模態(tài)包括兩種表達形式，分別由Maeda等[15,17]推導(dǎo)得到。對于一些難以求出自然模態(tài)的復(fù)雜結(jié)構(gòu)來說，Eatock Taylor等[18-19]推出了一種非常簡單的振型函數(shù)表達式：

(14)

式中：f1(x)表示剛體振型模態(tài)，本文采取此振型函數(shù)進行儲油船的水彈性分析。

根據(jù)振型函數(shù)的表達形式，可得儲油船的運動響應(yīng)函數(shù)w(x,y)表達式：

(15)

式中：ζj為j模態(tài)上的振型幅值。

考慮到線性勢流理論的假定，復(fù)速度勢φ包括入射波產(chǎn)生的入射勢、結(jié)構(gòu)存在形成的繞射勢以及物體運動產(chǎn)生的輻射勢之和，表示如下：

φ(x,y,z)=φI(x,y,z)+
φD(x,y,z)+φR(x,y,z)

(16)

(17)

(18)

2.2 速度勢及水動壓力分析

關(guān)于繞射勢和輻射勢邊界值問題的數(shù)值求解，通常采用格林函數(shù)法。將繞射勢φD和輻射勢φj(j=1,….N)分別與格林函數(shù)代入到格林第二公式中，得邊界積分方程：

(19)

其中(x0,y0,z0)和(x,y,z)分別為源點和場點坐標，G為滿足除物面條件外其他所有邊界條件的有限水深格林函數(shù)：

(20)

v=ω2/g為深水波數(shù)，

k0為色散方程ω2=gktanh(kh)的根，

J0為第一類貝塞爾函數(shù)。

積分方程(19)經(jīng)過物面單元離散，可表示為離散線性方程組形式：

[A]{φj}={B}

(21)

由于儲油船關(guān)于x-z和y-z兩平面對稱，此時系數(shù)矩陣及向量可表示為:

(22)

圖3 對稱區(qū)域簡化示意圖

各參數(shù)定義如圖3所示，對于矩陣A中各元素，滿足以下對稱關(guān)系：

A11=A22=A33=A44

(23)

A12=A21=A34=A43

(24)

A13=A31=A24=A42

(25)

A14=A41=A23=A32

(26)

為減小矩陣維數(shù)，可采用以下轉(zhuǎn)換:

(27)

其中矩陣R滿足以下關(guān)系：

(28)

這樣將式(27)和(28)代入線性方程組(21)可得：

(29)

根據(jù)矩陣[A]和[R]的特征，矩陣[A]是一個塊對角矩陣，則方程(29)為:

(30)

進一步將線性方程組(30)轉(zhuǎn)化為四個區(qū)域的求解形式，各區(qū)域表達式如下：

區(qū)域1：

(31)

(32)

(33)

(34)

區(qū)域2：

(35)

(36)

(37)

(38)

區(qū)域3：

(39)

(40)

(41)

(42)

區(qū)域4：

(43)

(44)

(45)

(46)

求得速度勢后，儲油船底板所受的壓力分布p(x,y)為：

p(x,y)=iωρ[φI+φD+φR]-ρgw(x,y)

(47)

2.3 儲油船的水彈性運動方程求解

由邊界元法求出速度勢及底板壓力分布后，將(15)和(47)式代入薄板的振動微分方程(10)可得：

(48)

對上述方程采用Galerkin積分法，即方程兩邊同時乘以fi(i=1,2,…,N)，并沿儲油船底部積分可得：

(49)

3 數(shù)值結(jié)果分析

根據(jù)上述數(shù)值方法，本文采用Fortran語言進行程序開發(fā)，并將計算結(jié)果與已知實驗和數(shù)值結(jié)果進行比較，驗證了本文程序的可行性。進一步利用所開發(fā)的程序，對儲油船在不同水深及波長下的水彈性響應(yīng)及變形、彎矩、廣義波浪激勵力和廣義水動力系數(shù)進行了分析。

3.1 程序的驗證

為了驗證本文數(shù)值方法及程序的可行性，選用Maeda等[15]實驗?zāi)Ｐ妥鳛樗憷⒂嬎憬Y(jié)果與實驗結(jié)果以及文獻[13]的數(shù)值結(jié)果進行比較，如圖4所示。實驗?zāi)Ｐ偷闹饕獏?shù)如表1所示。為了保證數(shù)值收斂和準確性，通常1個波浪長內(nèi)至少要包含14個單元，因此模型沿板長方向劃分60個單元，板寬方向16個單元，板高方向1個單元。

表1 實驗?zāi)Ｐ椭饕獏?shù)

圖4為波長分別2 m和4 m時，薄板左側(cè)即迎浪端的水彈性響應(yīng)位移最大時的結(jié)果。由圖可見，本文計算結(jié)果與實驗及文獻[5]結(jié)果基本吻合，只是在兩端位移數(shù)值上，本文計算結(jié)果稍大。且由圖可觀察到，波長較小時結(jié)構(gòu)變形程度要大于波長較大時的變形程度。

圖4 波長為2 m和 4 m時薄板水彈性位移分布

3.2 廣義波浪激勵力和廣義水動力系數(shù)的計算

本文所采用的儲油船數(shù)值模型按原型比尺1∶100取參數(shù)值：長L=4 m，寬B=0.9 m，吃水d=0.25 m，水深h=0.4 m、0.7 m、1 m，泊松比v=0.3，彎曲剛度D=50.832 Nm。型長方向劃分100個單元，型寬方向劃分20個單元，水深方向劃分5個單元。對不同水深條件下，運動方程(49)中儲油船所受的廣義波浪激勵力Fi隨相對波長變化進行了分析(取前6階模態(tài))，如圖5所示。由圖可知，隨著水深的增加，同階模態(tài)激勵力變化趨勢基本一致，且1、3階模態(tài)激勵力基本不變，其他階模態(tài)激勵力緩慢降低。

圖5 水深0.4 m,0.7 m和1.0 m時廣義波浪激勵力Fi計算結(jié)果圖

圖6 水深0.4 m時儲油船廣義附加質(zhì)量和輻射阻尼計算結(jié)果圖

此外，對不同水深條件下，運動方程(49)中儲油船的廣義附加質(zhì)量Maii和輻射阻尼Bbii隨波長的變化也進行了分析，如圖6～8所示。其中橫縱坐標取無量綱形式，廣義附加質(zhì)量和輻射阻尼都取到前6階模態(tài)。由圖可知，不同水深下的水動力系數(shù)變化趨勢也不相同：當(dāng)水深較淺時廣義附加質(zhì)量表現(xiàn)為低頻特征，廣義輻射阻尼表現(xiàn)為波頻特征；隨著水深的增加，廣義附加質(zhì)量也逐漸表現(xiàn)為波頻特征；當(dāng)水深增加到一定程度時，水動力系數(shù)表現(xiàn)為波頻和高頻兩種特征，因此水深的選取對儲油船的安全穩(wěn)定至關(guān)重要。從圖中可發(fā)現(xiàn)前6階模態(tài)水動力系數(shù)的累加基本可近似表達總的水動力系數(shù)值，且第1階模態(tài)即剛性模態(tài)水動力系數(shù)所占比例最大，第2階模態(tài)的水動力系數(shù)所占比例較小，第4階模態(tài)的廣義附加質(zhì)量在無量綱波數(shù)0.5～1.25時高于第3階模態(tài)值，其余模態(tài)水動力系數(shù)依模態(tài)階數(shù)降低。

3.3 不同水深下儲油船的水彈性響應(yīng)及變形

進一步研究不同水深條件下，儲油船在波長分別為0.8 m和1.0 m時的水彈性響應(yīng)，結(jié)果如圖9所示。由圖可得，當(dāng)水深從0.4 m～0.7 m變化時，儲油船水彈性響應(yīng)逐漸減小，但當(dāng)水深增大到1.0 m時水彈性響應(yīng)突然增大，并且水深為0.4 m～0.7 m時，結(jié)構(gòu)水彈性響應(yīng)對波長的敏感度小于水深1.0 m時結(jié)果。另外，由圖可知儲油船的迎浪端水彈性變形比較大，這一點在設(shè)計時必須引起注意。

本文計算了整個儲油船底板的彈性變形，下面給出了隨著水深變化，儲油船迎浪端出現(xiàn)最大響應(yīng)時的彈性變形圖，如圖10～13所示，其中左邊為波長0.8 m時計算結(jié)果，右邊是波長1.0 m計算結(jié)果。由圖可發(fā)現(xiàn)，水深為0.4～0.7 m時，結(jié)構(gòu)水彈性響應(yīng)對波長的敏感度要小于水深為1.0 m時的結(jié)果。

圖9 波長分別為0.8 m和1.0 m時儲油船的水彈性響應(yīng)

圖10 水深0.4 m時儲油船底板變形圖

圖11 水深0.55 m時儲油船底板變形圖

圖12 水深0.7 m時儲油船底板變形圖

圖13 水深1.0 m時儲油船底板變形圖

圖14 波長為0.8 m和1.0 m時儲油船的彎矩響應(yīng)分布

3.4 儲油船彎矩響應(yīng)分布

根據(jù)方程(13)計算不同水深時整個儲油船的彎矩分布，如圖14所示，其中上圖為波長0.8 m時的彎矩響應(yīng)分布，下圖為波長1.0 m時的彎矩響應(yīng)分布。由圖可知，當(dāng)水深從0.4～0.7 m變化時，整個儲油船的彎矩分布逐漸降低，到水深0.7 m時基本穩(wěn)定，但當(dāng)水深增加到1.0 m時，彎矩又突然增大，因此0.7 m水深可作為儲油船錨泊地點的相對理想水深；對于相同水深時，波長較大時的彎矩分布峰值大于波長較小的情況。

4 結(jié) 論

本文采用模態(tài)分析法中的模態(tài)展開法，將儲油船的水彈性響應(yīng)及輻射勢用一系列模態(tài)疊加表示，再應(yīng)用基于三維勢流理論的邊界積分方程和薄板振動微分方程求出各模態(tài)幅值，最后回代入儲油船的運動表達式得到水彈性。主要得出以下結(jié)論：

(1)通過與已知實驗及數(shù)值結(jié)果比較，驗證了本文方法能夠比較準確地預(yù)報超大型浮體水彈性響應(yīng)。

(2)隨著水深的增加，儲油船所受的同階模態(tài)廣義波浪激勵力變化趨勢基本一致，但廣義水動力系數(shù)變化趨勢并不相同。廣義附加質(zhì)量隨著水深的增加由低頻特征變?yōu)椴l特征最后同時具有波頻和高頻特征，廣義輻射阻尼由波頻特征變?yōu)橥瑫r具有波頻和高頻特征。

(3)隨著水深的增加，儲油船的水彈性響應(yīng)、整個結(jié)構(gòu)變形和彎矩響應(yīng)分布均隨之減小，到一定水深時趨于穩(wěn)定，但當(dāng)水深繼續(xù)增大時，水彈性響應(yīng)和彎矩會增加，因此超大型浮體系泊水深的選取對于水彈性響應(yīng)計算具有重要影響。

(4)整個儲油船的變形程度隨著入射波波長的增加而減小，這是由于船底部流體壓力分布的變化程度隨著波長的增加而減小造成的。

參 考 文 獻

[1]李文龍，劉亞東，譚家華. 基于模糊數(shù)學(xué)理論的我國戰(zhàn)略石油儲備方式?jīng)Q策[J].上海交通大學(xué)學(xué)報, 2004, 38(11)：1892-1896.

LI Wen-long, LIU Ya-dong, TAN Jia-hua. Decision-making for the facilities of strategic pectroleum reserve in china based on fuzzy mathematics theory[J]. Journal of Shanghai Jiao Tong University,2004,38(11):1892-1896.

[2]崔維成. 超大型海洋浮式結(jié)構(gòu)物水彈性響應(yīng)預(yù)報的研究現(xiàn)狀和發(fā)展方向[J]. 船舶力學(xué)，2002，6(1)：73-90.

CUI Wei-cheng. Current status and future direction in predicting the hydroelastic response of very large floating structures[J]. Journal of Ship Mechanics,2002,6(1):73-90.

[3]Bishop R E D, Price W G, Wu Y S. A general linear hydroelasticity theory of floating structures moving in a seaway. Phil[J]. Trans. R. Soc. Of Londan (series)A. 1986, 316:375-426.

[4]Wu C， Watanabe E， Utsunomiya T. An eigenfunction expansion-matching method for analyzing the wave-induced responses of an elastic floating plate[J]. Applied Ocean Research, 1995, 17(5):301-310.

[5]Hong D C, Hong S Y. Hydroelastic responses and drift forces of a very-long floating structure equipped with a pin-connected oscillating-water-column breakwater system[J]. Ocean Engineering, 2007，34(5-6):696-708.

[6]Murai M, Kagemoto H, Fujino M. On the hydroelastic responses of a very large floating structure in waves[J]. Marine Science and Technology, 1999, 4(3):123-153.

[7]Fu S X, Moan T, Chen X J, et al. Hydroelastic analysis of flexible floating interconnected structures[J]. Ocean Engineering, 2007, 34(11-12):1516-1531.

[8]Watanabe E, Utsunomiya T, Wang C M, et al. Benchmark hydroelastic responses of a circular VLFS under wave action[J]. Engineering Structures, 2006, 28(3):423-430.

[9]Pham D C, Wang C M, Utsunomiya T. Hydroelastic analysis of pontoon-type circular VLFS with an attached submerged plate[J]. Applied ocean Research, 2008, 30(4):287-296.

[10]Wu Y S. Hydroelasticity of floating bodies[D]. Brunel University, UK,1984.

[11]金晶哲，崔維成，劉應(yīng)中. 預(yù)報超大型浮體水彈性響應(yīng)的模態(tài)函數(shù)展開方法和特征函數(shù)展開方法比較[J]. 船舶力學(xué)，2003，7(4)：86-98.

JIN Jin-zhe, CUI Wei-cheng, LIU Ying-zhong. Comparison of modal function expansion method with eigenfunction expansion method for prediction of hydroelastic responses of VLFS[J]. Journal of Ship Mechanics,2003,7(4):86-98.

[12]閆紅梅，崔維成，劉應(yīng)中.超大型浮體水彈性分析的平板林函數(shù)法[J].海洋工程，2003, 21(4)：8-21.

YAN Hong-mei, CUI Wei-cheng, LIU Ying-zhong. Hydroelastic analysis of VLFS using Plate Green Function Method[J]. The Ocean Engineering, 2003, 21(4):8-21.

[13]滕 斌，勾 瑩.大型浮體水彈性作用的頻域分析[J].工程學(xué)，2006, 23(s2)：36-47.

TENG Bin, GOU Ying. Hydroelastic analysis of very large floating structure in frequency domain[J]. Engineering mechanics,2006,23(增2):36-47.

[14]張淑華，韓滿生.用直接法分析超大型浮體的水彈性響應(yīng)[J].海洋工程，2004, 22(1)：9-18.

ZHANG Shu-hua, HAN Man-sheng. Analysis of hydroelastic response of VLFS by use of the direct method[J].The Ocean Engineering, 2004, 22(1):9-18.

[15]Maeda H, Masuda K, Miyajima S,et al. Hydrielastic reponsesof pontoon type very large floating offshore structure[J]. Jour of Soc Nav Archit Japan, 1995, 178：203-212.

[16]李文龍，譚家華.我國戰(zhàn)略石油儲備船系泊系統(tǒng)安全性評估[J].中國造船，2005,46(2):87-94.

LI Wen-long, TAN Jia-hua. Safety assessment of mooring system of oil storage vessels for strategic petroleum reserve in china [J].Shipbuding of China,2005, 46(2):87-94.

[17]Gran S. A course in ocean engineering[M]. Elsevier, Amsterdam,1992.

[18]Eatock Taylor R, Ohkusu M. Green functions for hydroelasticanalysis of vibrating free-free beams and plates[J]. Applied Ocean Research, 2000, 22：295-314.

[19]Eatock Taylor R. Wet or dry modes in linear hydroelasticity-why modes [C]. Hydroelasticity in Marine Technology, Oxford, UK, 2003：239-250.