金 鑫，劉飛飛，鄒寅劼，任 律，劉名名,6，李登松,朱紅鈞

(1.成都理工大學(xué) 能源學(xué)院，成都 610059； 2.河海大學(xué) 海岸災(zāi)害及防護(hù)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210098；3.天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072；4.淄博市公用事業(yè)服務(wù)中心，淄博 255000；5.水電水利規(guī)劃設(shè)計(jì)總院，北京 100120； 6.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連 116024； 7.四川農(nóng)業(yè)大學(xué) 水利水電學(xué)院，雅安 625014；8.西南石油大學(xué) 油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室，成都 610500)

經(jīng)濟(jì)的發(fā)展離不開能源的支撐，隨著陸上油氣資源的日益枯竭，各國均加大海洋油氣資源的開采和勘察力度。船載運(yùn)輸具有成本低的優(yōu)點(diǎn)，目前被廣泛采用，流體晃蕩問題也不可避免，同時(shí)也是海洋工程中的熱點(diǎn)問題。在一些極端海況下，船載流體在多自由度海浪作用下極易不穩(wěn)定，形成的劇烈晃蕩波極有可能破壞液艙的結(jié)構(gòu)甚至引起運(yùn)輸船的傾覆，因此亟需對(duì)流體晃蕩問題深入研究。

早期研究以勢(shì)流理論為基礎(chǔ)，對(duì)單自由度縱蕩、升沉、縱蕩與升沉耦合激勵(lì)下的流體晃蕩問題進(jìn)行分析(Faltinsen[1]；Benjamin和Ursell[2]；Hill[3]；Frandsen[4])，用于預(yù)測(cè)晃蕩波高及頻譜特性。Faltinsen等[5]證實(shí)高階模態(tài)對(duì)晃蕩波的非線性波浪特性具有重要影響。Frandsen和Peng[6]試驗(yàn)地研究單獨(dú)縱蕩或升沉激勵(lì)下的流體晃蕩問題，發(fā)現(xiàn)充液深度對(duì)波浪的非線性響應(yīng)具有重要貢獻(xiàn)。衛(wèi)志軍等[7]對(duì)超大型儲(chǔ)船內(nèi)的流體晃蕩進(jìn)行實(shí)驗(yàn)研究，針對(duì)晃蕩荷載進(jìn)行系統(tǒng)的分析。Xue和Lin[8]和Xue等[9]對(duì)多種激勵(lì)形式下的晃蕩波進(jìn)行分析并討論內(nèi)部結(jié)構(gòu)對(duì)共振波的影響機(jī)制。

隨著計(jì)算機(jī)硬件的發(fā)展，數(shù)值模擬逐步流行。早期的數(shù)值模型以Euler方程為基礎(chǔ)，相關(guān)學(xué)者對(duì)二維和三維問題進(jìn)行系統(tǒng)的研究，數(shù)值方法包括邊界元法、有限差分法、光滑粒子流體動(dòng)力學(xué)法、移動(dòng)粒子半隱式法等(Nakayama和Washizu[10]；Okamoto和Kawahara[11]；Cho等[12]；Gingold和Monaghan[13]；Yang等[14]；Wu等[15]；田昀艷[16])。由于勢(shì)流模型無法真實(shí)反映完整的流體運(yùn)動(dòng)規(guī)律，Navier-Stokes模型逐步發(fā)展起來，研究人員對(duì)強(qiáng)非線性和完全非線性晃蕩問題進(jìn)行數(shù)值模擬研究，對(duì)晃蕩波高與激勵(lì)頻率的響應(yīng)關(guān)系進(jìn)行細(xì)致分析，討論激勵(lì)頻率、激勵(lì)振幅、充液深度對(duì)晃蕩波高與荷載的影響規(guī)律(Ramaswamy[17]；孫江龍和葉恒奎[18]；Shao等[19]；Luo等[20-21])。

Frandsen[4]推導(dǎo)出縱蕩與升沉耦合激勵(lì)下流體晃蕩問題的解析解，發(fā)現(xiàn)耦合激勵(lì)頻率的線性組合滿足一定條件時(shí)，波高會(huì)出現(xiàn)指數(shù)型增長(zhǎng)。Ning等[22]開發(fā)了一套基于邊界元算法的數(shù)值模型證實(shí)Frandsen[4]的結(jié)論。Sriram等[23]通過有限元模型分析耦合激勵(lì)下的流體晃蕩問題，得到流體晃蕩響應(yīng)的頻譜特性。相關(guān)研究證實(shí)共振波有激發(fā)的趨勢(shì)，但長(zhǎng)時(shí)間的試驗(yàn)或數(shù)值模擬結(jié)果很少被研究，共振波的具體形式或頻率組合規(guī)律的研究相對(duì)較少。本文將著重分析耦合激勵(lì)下共振波及其非線性特性，同時(shí)也將對(duì)耦合共振激勵(lì)的頻率組合規(guī)律開展研究工作。

1 線性晃蕩理論

考慮長(zhǎng)度為L(zhǎng)、水深為h的二維棱柱形封閉充液箱在水平x方向(縱蕩)和垂直y方向(升沉)激勵(lì)下的受迫運(yùn)動(dòng)，箱體運(yùn)動(dòng)軌跡分別為

X=X(t)=ahcos(ωht)
Z=Z(t)=avcos(ωvt)

(1)

式中：ah和av分別為水平與垂向的激勵(lì)振幅(m)，ωh和ωv分別代表水平與垂向的激勵(lì)頻率，rad/s。

假設(shè)流體是無粘、不可壓，則流體是有勢(shì)φ，其速度勢(shì)φ與波高ζ可以表示為

(2)

式中：kn=(2n+1)π/L。

經(jīng)過數(shù)學(xué)推導(dǎo)可得(具體過程參考Frandsen[4])

(3)

因此，得到縱蕩與升沉耦合激勵(lì)下流體晃蕩的一般形式

(4)

求解上述非齊次方程，疊加不同模態(tài)下的波高，即可得到縱蕩與升沉耦合激勵(lì)下晃蕩問題的理論波高。如果只考慮升沉激勵(lì)，不考慮縱蕩的影響，可得到齊次Mathieu方程，即為單獨(dú)升沉激勵(lì)晃蕩問題。當(dāng)該方程的解呈現(xiàn)出指數(shù)型增長(zhǎng)，即升沉參數(shù)位于不穩(wěn)定區(qū)，其他情況下流體的運(yùn)動(dòng)強(qiáng)度微弱，屬于穩(wěn)定區(qū)內(nèi)的流體晃蕩。根據(jù)Frandsen[4]，第一與第二不穩(wěn)定區(qū)分別位于ω1/ωn= 0.5與ω1/ωn= 1.0的附近。

2 數(shù)值模型與驗(yàn)證

2.1 流體運(yùn)動(dòng)控制方程

本文采用不可壓縮流體運(yùn)動(dòng)控制方程為Navier-Stokes方程組，可表示為

(5)

(6)

式中：i和j為自由指標(biāo)(i,j= 1, 2, 3)，分別代表x,y,z三個(gè)方向上的分量；ui為流體速度分量；ρ為流體密度；p為壓強(qiáng)；gi為i方向的重力加速度分量；fi為i方向的箱體運(yùn)動(dòng)加速度分量。由于本文只考慮縱蕩和升沉運(yùn)動(dòng)，因此流體運(yùn)動(dòng)控制方程中fi可寫為

(7)

2.2 數(shù)值模型驗(yàn)證

Frandsen[4]利用有限差分?jǐn)?shù)值模型研究縱蕩與升沉耦合激勵(lì)下的晃蕩問題。箱體長(zhǎng)度L=1.0 m，水深h=0.5 m，一階自然頻率ω1=5.316 6 rad/s?？v蕩與升沉均按余弦形式做周期往復(fù)運(yùn)動(dòng)，參數(shù)分別為：ah= 0.000 506 m，ωh= 0.98ω1，av= 0.272 4 m和ωh= 0.8ω1。二維計(jì)算域采用400×1×240的均勻網(wǎng)格系統(tǒng)，最小網(wǎng)格為0.002 5 m。圖 1給出左邊壁自由液面高程的時(shí)間歷程的理論解、文獻(xiàn)數(shù)據(jù)與本模型結(jié)果的比較關(guān)系及相應(yīng)的波高相位圖。從圖1-a中可以發(fā)現(xiàn)：本模型與Frandsen的數(shù)值解[4]吻合良好，二者與理論解在無量綱時(shí)間t×ω1= 52之后偏差逐步增大，主要由于理論解忽略波浪的非線性?？紤]到縱蕩激勵(lì)接近一階自振頻率，一階模態(tài)共振波被激發(fā)。理論解與本模型結(jié)果的波高相位圖分別如圖1-b和圖1-c所示，可以發(fā)現(xiàn)，理論解呈現(xiàn)出對(duì)稱螺旋狀的增長(zhǎng)，而本模型考慮流體運(yùn)動(dòng)的非線性，表現(xiàn)出明顯的不對(duì)稱性。

1-a 時(shí)程曲線對(duì)比

1-b 理論解軌跡1-c 數(shù)值解軌跡注:ah = 0.000 506 m,ωh = 0.98 ω1, av = 0.272 4 m,ωh = 0.8 ω1。圖1 縱蕩與升沉耦合激勵(lì)下箱體左邊壁自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖Fig.1 Time series of free surface displacements at the left wall and corresponding wave-phase diagrams under coupled surge and heave excitations

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

在實(shí)際情況中，運(yùn)輸液艙在單獨(dú)縱蕩或升沉激勵(lì)作用下的情況較少，一般以耦合激勵(lì)為主。本文主要討論縱蕩與升沉耦合激勵(lì)下的流體晃蕩問題，分為三類，Case 1和Case 2：縱蕩激勵(lì)頻率低于一階自然頻率且升沉處于穩(wěn)定區(qū)；Case 3：縱蕩激勵(lì)頻率高于一階自然頻率且升沉處于穩(wěn)定區(qū)；Case 4：縱蕩激勵(lì)頻率遠(yuǎn)離一階自然頻率，升沉參數(shù)處于第一不穩(wěn)定區(qū)。相關(guān)參數(shù)總結(jié)如表1。鑒于Case 4必然會(huì)激發(fā)共振晃蕩波，下面重點(diǎn)討論Case 1～Case 3，在此基礎(chǔ)上將討論耦合激勵(lì)晃蕩自由面高程與激勵(lì)頻率間的響應(yīng)關(guān)系。箱體尺寸為：長(zhǎng)度L= 0.5 m，水深h= 0.3 m，一階自然頻率ω1= 7.672 rad/s。

表1 縱蕩與升沉耦合激勵(lì)參數(shù)Tab.1 Excitation parameters for coupled surge and heave sloshing

3.1 頻率和值為ω1

Case 1和Case 2的縱蕩頻率偏離一階自然頻率，均屬于非共振激勵(lì)；升沉激勵(lì)參數(shù)位于穩(wěn)定區(qū)，二者均無法單獨(dú)激發(fā)明顯的波浪運(yùn)動(dòng)。激勵(lì)頻率的和值ωh+ωv等于一階自然頻率ω1，屬于共振頻率。Case 1左邊壁處自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖如圖2所示。盡管有指數(shù)增長(zhǎng)的趨勢(shì)，由于縱蕩頻率遠(yuǎn)低于共振頻率且激勵(lì)振幅較小，共振特性不明顯。理論與數(shù)值結(jié)果均呈現(xiàn)出增長(zhǎng)包絡(luò)線，如圖 2-a所示，這與單獨(dú)縱蕩激勵(lì)晃蕩較為接近；隨著波高的逐步增大，數(shù)值解逐步表現(xiàn)出非對(duì)稱性，如波高相位圖如圖 2-c所示。

2-a 時(shí)程曲線對(duì)比

2-b 理論解軌跡2-c 數(shù)值解軌跡注:ah=0.02 m,ωh=0.2 ω1, av=0.130 2 m,ωv=0.8 ω1。圖2 Case 1左邊壁自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖Fig.2 Time series of free surface displacements at the left wall of Case 1 and corresponding wave-phase diagrams under coupled surge and heave excitations

經(jīng)過長(zhǎng)時(shí)間模擬Case 2，得到左邊壁處自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖，如圖 3所示?？梢园l(fā)現(xiàn)，由于非線性的影響，波高不會(huì)無限增長(zhǎng)，且包絡(luò)周期明顯減小。共振特征十分明顯，這與單獨(dú)縱蕩或穩(wěn)定升沉激勵(lì)下的晃蕩有著顯著不同。Case 2與Case 1比較而言，縱蕩頻率逐步靠近共振頻率(一階自然頻率)，且縱蕩加速度逐步增大，波高增長(zhǎng)明顯，且能達(dá)到縱蕩振幅的3.5倍，接近7.0 cm，非線性逐步占據(jù)主導(dǎo)，因此理論解與數(shù)值解在無量綱時(shí)間t×ω1= 250之后的偏差逐步增大，數(shù)值解的不對(duì)稱性十分明顯，如圖3-c的波高相位圖所示，而理論解的波高相位圖(圖3-b)始終呈現(xiàn)出對(duì)稱性。不同時(shí)刻自由液面波形如圖 4所示，可以發(fā)現(xiàn)，明顯的共振晃蕩波被激發(fā)，最大波高總出現(xiàn)在箱體壁面，波長(zhǎng)接近箱體長(zhǎng)度的2倍，明顯的一階共振波形。

3-a 時(shí)程曲線對(duì)比

3-b 理論解軌跡3-c 數(shù)值解軌跡注:ah = 0.02 m,ωh = 0.45 ω1, av = 0.040 3 m,ωv = 0.55 ω1。圖3 Case 2左邊壁自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖Fig.3 Time series of free surface displacements at the left wall of Case 2 and corresponding wave-phase diagrams under coupled surge and heave excitations

圖4 Case 2自由面波形圖Fig.4 Free surface snapshots of Case 2

3.2 頻率差值為ω1

激勵(lì)頻率差值ωv-ωh等于一階自然頻率ω1的Case 3也進(jìn)行了數(shù)值模擬研究，相應(yīng)左邊壁處自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖如圖 5所示。與Case 2比較，Case 3的縱蕩頻率較低，縱蕩振幅一致，盡管升沉激勵(lì)均處于穩(wěn)定區(qū)，但Case 3的波高較Case 2大，表明差值組合較和值組合的工況更容易激發(fā)共振波形。同樣地，由于非線性的影響，數(shù)值波高與線性解在波高較大時(shí)差異明顯，呈現(xiàn)出非對(duì)稱的形狀且大于理論解，如圖 5所示；非線性特性也可以從自由面高程相位圖(圖 5-c)發(fā)現(xiàn)。不同時(shí)刻自由液面波形如圖6所示，可以發(fā)現(xiàn)，明顯的一階模態(tài)的共振晃蕩波被激發(fā)，與Case 2類似，波長(zhǎng)接近箱體長(zhǎng)度的2倍。

5-a 時(shí)程曲線對(duì)比

5-b 理論解軌跡5-c 數(shù)值解軌跡注:ah = 0.02 m,ωh = 0.353 3 ω1, av = 0.037 5 m,ωv = 1.333 3 ω1。圖5 Case 3左邊壁自由面高程的時(shí)間歷程曲線及相應(yīng)的波高相位圖Fig.5 Time series of free surface displacements at the left wall of Case 3 and corresponding wave-phase diagrams under coupled surge and heave excitations

圖6 Case 3自由面波形圖Fig.6 Free surface snapshots of Case 3

從Case 1～Case 3的模擬結(jié)果可以發(fā)現(xiàn)，即使單自由度的穩(wěn)定區(qū)升沉激勵(lì)或縱蕩非共振激勵(lì)不能激發(fā)出連續(xù)增長(zhǎng)的共振形態(tài)，當(dāng)耦合激勵(lì)滿足一定條件時(shí)，即：它們的差值或和值等于自然頻率時(shí)，也能激發(fā)共振晃蕩波。

3.3 掃頻響應(yīng)

上述研究表明：當(dāng)縱蕩和升沉激勵(lì)均處于非共振時(shí)，耦合激勵(lì)下的晃蕩在某些組合下也會(huì)出現(xiàn)共振特性，下面將對(duì)縱蕩和升沉耦合激勵(lì)下的流體晃蕩響應(yīng)進(jìn)行系統(tǒng)的分析。在保證縱蕩和升沉激勵(lì)加速度恒定的情況下，系統(tǒng)地改變縱蕩和升沉激勵(lì)頻率，考察波高η隨頻率的響應(yīng)過程，即：掃頻分析，分為兩種工況：第一種，縱蕩激勵(lì)非共振，系統(tǒng)地改變升沉激勵(lì)參數(shù)；第二種，升沉激勵(lì)處于穩(wěn)定區(qū)，系統(tǒng)地改變縱蕩激勵(lì)參數(shù)。

3.3.1 縱蕩掃頻響應(yīng)

升沉運(yùn)動(dòng)在數(shù)值模型中設(shè)置為av=0.130 2 m，ωv= 0.8ω1，屬于穩(wěn)定區(qū)的升沉晃蕩。保持升沉參數(shù)不變，系統(tǒng)地改變縱蕩激勵(lì)頻率，相應(yīng)的角頻率為0.383 6 rad/s到17.645 6 rad/s，對(duì)應(yīng)[0.05ω1, 2.3ω1]。保持縱蕩加速度不變，即：ahωh2/g=0.01 g的情況下，水平激勵(lì)頻率從0.05ω1系統(tǒng)地增大2.25ω1，對(duì)應(yīng)縱蕩激勵(lì)振幅為[0.666 7 m, 0.000 3 m]。理論解的頻率間隔為0.01ω1，數(shù)值解的頻率間隔為0.05ω1。圖7給出自由面高程-縱蕩激勵(lì)頻率的響應(yīng)曲的理論解與數(shù)值模擬結(jié)果的對(duì)比關(guān)系?？梢园l(fā)現(xiàn)，當(dāng)激勵(lì)頻率線性組合偏離自然頻率時(shí)，理論解和數(shù)值結(jié)果較為吻合，當(dāng)頻率線性組合接近自然頻率時(shí)，數(shù)值結(jié)果和解析解存在一定的差異，主要是由于共振誘發(fā)的非線性引起的。圖7中各峰值對(duì)應(yīng)頻率分別為0.18ω1、0.625ω1、1.0ω1、1.45ω1、1.75ω1和2.25ω1，其中：0.18ω1和升沉頻率的和值接近一階自然頻率ω1；0.625ω1和升沉頻率的和值接近第二階自然頻率ω2；一階自然頻率ω1，其自由面高程最大；1.45ω1接近第二階自然頻率ω2；1.75ω1接近第三階自然頻率接近ω3；2.25ω1和升沉激勵(lì)頻率的差值接近第二階自然頻率ω2。結(jié)果表明：當(dāng)二者頻率和值或差值接近自然頻率時(shí)，自由面高程明顯高于低于該值的區(qū)域，單獨(dú)縱蕩激勵(lì)接近自然頻率時(shí)，也會(huì)激發(fā)共振波。

圖7 穩(wěn)定區(qū)升沉激勵(lì)與掃頻縱蕩下左邊壁自由面高程-縱蕩激勵(lì)頻率相應(yīng)曲線Fig.7 Free surface displacement vs. surge frequency for the sweeping surge test with stable surge frequency

3.3.2 升沉掃頻響應(yīng)

縱蕩參數(shù)固定不變，設(shè)置為：ah= 0.005 m，ωh= 0.5ω1。保持升沉加速度不變，即：avωv2/g=0.25g，升沉頻率從0.05ω1到 2.3ω1系統(tǒng)地改變，對(duì)應(yīng)激勵(lì)振幅從16.666 8 m遞減到0.007 9 m。理論解的頻率間隔為0.01ω1，數(shù)值解的頻率間隔為0.05ω1。圖8給出左邊壁自由面高程-升沉激勵(lì)頻率的響應(yīng)曲線結(jié)果，同時(shí)也標(biāo)出升沉激勵(lì)所在不穩(wěn)定區(qū)的分布。鑒于線性理論解在不穩(wěn)定區(qū)會(huì)指數(shù)型增長(zhǎng)，無量綱自由面高程非常大，不便于結(jié)果對(duì)比，因此在圖中單獨(dú)列出比較完整的解析響應(yīng)(圖8)，其中虛線覆蓋區(qū)表示自由面高程迅速增長(zhǎng)的區(qū)域?？梢园l(fā)現(xiàn)在第一不穩(wěn)定區(qū)自由面高程明顯大于其他工況，同時(shí)也可以看出，橫坐標(biāo)為0.5、1.0、1.5、1.75、2和2.25時(shí)，晃蕩的波高明顯大于偏離該頻率的工況。對(duì)于橫坐標(biāo)為0.5工況，它和縱蕩激勵(lì)頻率的和等于一階自然頻率ω1；橫坐標(biāo)為1.0即是一階自然頻率ω1，同時(shí)也處于第二不穩(wěn)定區(qū)；橫坐標(biāo)為1.5時(shí)，它和縱蕩激勵(lì)頻率的差值等于一階自然頻率ω1；橫坐標(biāo)為1.75時(shí)，縱蕩激勵(lì)頻率接近第三階自然頻率ω3；橫坐標(biāo)為2.0時(shí)，處于第一不穩(wěn)定區(qū)，對(duì)應(yīng)的波高最大；橫坐標(biāo)為2.25時(shí)，它和升沉激勵(lì)頻率的差值接近第三階自然頻率ω3。

圖8 非共振縱蕩激勵(lì)與掃頻升沉下左邊壁自由面高程-升沉激勵(lì)頻率相應(yīng)曲線Fig.8 Free surface displacement vs. heave frequency for the sweeping heave test with off-resonant surge frequency

4 結(jié)論

本文通過有限差分法求解不可壓縮Navier-Stokes方程，建立縱蕩與升沉激勵(lì)下流體晃蕩的三維數(shù)值模型，數(shù)值模擬結(jié)果與理論解、文獻(xiàn)中的數(shù)據(jù)比較分析，證明本模型在模擬耦合激勵(lì)下晃蕩問題的準(zhǔn)確性。在此基礎(chǔ)上，通過系統(tǒng)的數(shù)值模擬研究，得出如下結(jié)論：

(1)耦合激勵(lì)頻率的差值或和值接近自然頻率時(shí)，共振晃蕩波可以被激發(fā)。

(2)升沉處于穩(wěn)定區(qū)的耦合晃蕩以縱蕩激勵(lì)晃蕩為主，當(dāng)二者激勵(lì)頻率和值或差值接近各階自然頻率時(shí)，會(huì)誘發(fā)共振；單獨(dú)縱蕩激勵(lì)接近自然頻率時(shí)，也會(huì)激發(fā)共振晃蕩波。

(3)非共振縱蕩的耦合激勵(lì)晃蕩以升沉激勵(lì)晃蕩為主，尤其是升沉位于不穩(wěn)定區(qū)的工況。