縱蕩與升沉耦合激勵下晃蕩波特性及響應規律研究

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

(1.成都理工大學 能源學院，成都 610059； 2.河海大學 海岸災害及防護教育部重點實驗室，南京 210098；3.天津大學 水利工程仿真與安全國家重點實驗室，天津 300072；4.淄博市公用事業服務中心，淄博 255000；5.水電水利規劃設計總院，北京 100120； 6.大連理工大學 海岸和近海工程國家重點實驗室，大連 116024； 7.四川農業大學 水利水電學院，雅安 625014；8.西南石油大學 油氣藏地質及開發工程國家重點實驗室，成都 610500)

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

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

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

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

1 線性晃蕩理論

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

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

(1)

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

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

(2)

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

經過數學推導可得(具體過程參考Frandsen[4])

(3)

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

(4)

求解上述非齊次方程，疊加不同模態下的波高，即可得到縱蕩與升沉耦合激勵下晃蕩問題的理論波高。如果只考慮升沉激勵，不考慮縱蕩的影響，可得到齊次Mathieu方程，即為單獨升沉激勵晃蕩問題。當該方程的解呈現出指數型增長，即升沉參數位于不穩定區，其他情況下流體的運動強度微弱，屬于穩定區內的流體晃蕩。根據Frandsen[4]，第一與第二不穩定區分別位于ω1/ωn= 0.5與ω1/ωn= 1.0的附近。

2 數值模型與驗證

2.1 流體運動控制方程

本文采用不可壓縮流體運動控制方程為Navier-Stokes方程組，可表示為

(5)

(6)

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

(7)

2.2 數值模型驗證

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

1-a 時程曲線對比

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

3 計算結果與分析

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

表1 縱蕩與升沉耦合激勵參數Tab.1 Excitation parameters for coupled surge and heave sloshing

3.1 頻率和值為ω1

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

2-a 時程曲線對比

2-b 理論解軌跡2-c 數值解軌跡注:ah=0.02 m,ωh=0.2 ω1, av=0.130 2 m,ωv=0.8 ω1。圖2 Case 1左邊壁自由面高程的時間歷程曲線及相應的波高相位圖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

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

3-a 時程曲線對比

3-b 理論解軌跡3-c 數值解軌跡注:ah = 0.02 m,ωh = 0.45 ω1, av = 0.040 3 m,ωv = 0.55 ω1。圖3 Case 2左邊壁自由面高程的時間歷程曲線及相應的波高相位圖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

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

5-a 時程曲線對比

5-b 理論解軌跡5-c 數值解軌跡注:ah = 0.02 m,ωh = 0.353 3 ω1, av = 0.037 5 m,ωv = 1.333 3 ω1。圖5 Case 3左邊壁自由面高程的時間歷程曲線及相應的波高相位圖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的模擬結果可以發現，即使單自由度的穩定區升沉激勵或縱蕩非共振激勵不能激發出連續增長的共振形態，當耦合激勵滿足一定條件時，即：它們的差值或和值等于自然頻率時，也能激發共振晃蕩波。

3.3 掃頻響應

上述研究表明：當縱蕩和升沉激勵均處于非共振時，耦合激勵下的晃蕩在某些組合下也會出現共振特性，下面將對縱蕩和升沉耦合激勵下的流體晃蕩響應進行系統的分析。在保證縱蕩和升沉激勵加速度恒定的情況下，系統地改變縱蕩和升沉激勵頻率，考察波高η隨頻率的響應過程，即：掃頻分析，分為兩種工況：第一種，縱蕩激勵非共振，系統地改變升沉激勵參數；第二種，升沉激勵處于穩定區，系統地改變縱蕩激勵參數。

3.3.1 縱蕩掃頻響應

升沉運動在數值模型中設置為av=0.130 2 m，ωv= 0.8ω1，屬于穩定區的升沉晃蕩。保持升沉參數不變，系統地改變縱蕩激勵頻率，相應的角頻率為0.383 6 rad/s到17.645 6 rad/s，對應[0.05ω1, 2.3ω1]。保持縱蕩加速度不變，即：ahωh2/g=0.01 g的情況下，水平激勵頻率從0.05ω1系統地增大2.25ω1，對應縱蕩激勵振幅為[0.666 7 m, 0.000 3 m]。理論解的頻率間隔為0.01ω1，數值解的頻率間隔為0.05ω1。圖7給出自由面高程-縱蕩激勵頻率的響應曲的理論解與數值模擬結果的對比關系。可以發現，當激勵頻率線性組合偏離自然頻率時，理論解和數值結果較為吻合，當頻率線性組合接近自然頻率時，數值結果和解析解存在一定的差異，主要是由于共振誘發的非線性引起的。圖7中各峰值對應頻率分別為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和升沉激勵頻率的差值接近第二階自然頻率ω2。結果表明：當二者頻率和值或差值接近自然頻率時，自由面高程明顯高于低于該值的區域，單獨縱蕩激勵接近自然頻率時，也會激發共振波。

圖7 穩定區升沉激勵與掃頻縱蕩下左邊壁自由面高程-縱蕩激勵頻率相應曲線Fig.7 Free surface displacement vs. surge frequency for the sweeping surge test with stable surge frequency

3.3.2 升沉掃頻響應

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

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

4 結論

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

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

(2)升沉處于穩定區的耦合晃蕩以縱蕩激勵晃蕩為主，當二者激勵頻率和值或差值接近各階自然頻率時，會誘發共振；單獨縱蕩激勵接近自然頻率時，也會激發共振晃蕩波。

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