廖 斌，陳善群

(安徽工程大學 建筑工程學院，安徽 蕪湖 241000)

浮體一般是指漂浮于液體表面的物體，諸如浮標、船只、浮橋以及大型網(wǎng)箱等，在海洋及河流等大型水體中均極為常見。浮體在水體中易受波浪的驅(qū)動而產(chǎn)生一定的搖晃，當波浪作用較強時甚至會出現(xiàn)浮體的傾覆以及相互拍擊等現(xiàn)象，直接影響浮體的穩(wěn)定性，進一步危及浮體上設備的正常工作以及人員的安全[1]。同時，波浪中浮體運動的研究涉及到流固耦合、自由面精細追蹤等一系列學術(shù)問題。因此，研究波浪中的浮體運動具有較為重要的現(xiàn)實意義和學術(shù)價值。

國內(nèi)外對于波浪中浮體運動的研究由來已久，使用的手段目前主要集中在實驗和數(shù)值研究兩方面。Zhao[2]等對三種類型波浪(規(guī)則波、聚焦波和規(guī)則波與聚焦波的組合)驅(qū)動下倒T型箱型浮體運動進行了實驗觀測，獲取了浮體在3種類型波浪作用下的起伏位移以及搖蕩角度等相關(guān)數(shù)據(jù)。He[3]等使用電荷耦合器件攝像機捕獲浮體的實時運動，對6種波高和6種周期條件下波浪驅(qū)動中的箱型浮體運動的非線性行為進行了實驗研究。Wang[4]等對規(guī)則和不規(guī)則波浪驅(qū)動下沙漏型浮體運動進行了實驗研究，獲得了波浪作用下浮體運動的搖蕩周期、粘性阻尼和運動響應等實驗數(shù)據(jù)。Hashiura[5]等和Sun[6]等對長寬比較大的橋式浮體在規(guī)則波浪中的運動響應分別進行了實驗研究，重點對浮體所受的波浪力、彎矩等數(shù)據(jù)進行了深入分析。Ji[7]等對并列水平圓柱和矩形單浮橋浮式防波堤的水動力性能進行了實驗研究，對浮體的系泊索力和三自由度運動(搖擺，升沉和滾動)進行了系統(tǒng)評估。Yang[8]等設計了一種新型的壓載水浮式防波堤，并對該浮體在二維規(guī)則波驅(qū)動下的流體力學效率和運動響應進行了深入分析。顯然，實驗研究可以獲取浮體運動的大量直觀數(shù)據(jù)，但所需的測試場地以及采集設備的經(jīng)濟成本較高，能實現(xiàn)的波浪類型和浮體種類也相對有限。基于此，部分學者將目光聚焦于波浪驅(qū)動下浮體運動的數(shù)值研究，其中采用的數(shù)值方法主要分為有網(wǎng)格方法和無網(wǎng)格方法兩大類。采用有網(wǎng)格方法，Wu[9]等基于二維三叉樹方法生成運動網(wǎng)格，數(shù)值研究了具有漂浮結(jié)構(gòu)的拖曳罐與非線性波浪的相互作用。Hu[10]等研發(fā)了一套流體動力學(CFD)程序AMAZON-SC，并數(shù)值研究了規(guī)則波驅(qū)動下波能轉(zhuǎn)換器的運動響應規(guī)律。Jeong[11]等發(fā)展了一種基于改進標記密度法的數(shù)值模型來模擬浮體和自由表面的運動，計算的固定圓柱陣列與波浪的相互作用以及海上平臺的響應幅度與研究數(shù)據(jù)相符。鑒于有網(wǎng)格方法在處理流固耦合問題時網(wǎng)格易在流體自由面與固體交界面處產(chǎn)生畸變，需對網(wǎng)格進行調(diào)整或者重構(gòu)，耗費計算資源且影響計算精度。相較之下，無網(wǎng)格方法將固體和流體均離散成粒子，粒子之間的拓撲關(guān)系可不受物性參數(shù)的影響。因此，由粒子構(gòu)成的物質(zhì)界面隨時間更新時可避免網(wǎng)格的調(diào)整或者重構(gòu)，在提高計算效率的同時還可消除物質(zhì)傳遞帶來的不利影響。光滑粒子流體動力學方法(Smoothed Particle Hydrodynamics，SPH)以及移動粒子半隱式方法(Moving Particle Semi-implicit，MPS)為目前較為流行的兩種無網(wǎng)格方法。Omidvar[12]等使用具有嵌入式黎曼求解器的任意拉格朗日-歐拉公式實現(xiàn)可變粒子質(zhì)量分布，開發(fā)了一種可變質(zhì)量分布的SPH方法，并數(shù)值研究了二維楔形和三維圓錐形自由浮體隨波浪的運動規(guī)律。曹文瑾[13]等在MPS中建立了剛體與流體之間的單向及雙向耦合模型，并對活塞式造波池中長板型浮體運動進行了數(shù)值分析。

研究擬采用SPH方法建立一種適用于研究活塞式造波驅(qū)動下浮體運動的數(shù)值模型，并獲取箱式浮體六自由度條件下運動數(shù)據(jù)的變化曲線，發(fā)掘波浪驅(qū)動下箱式浮體運動的響應規(guī)律，為浮體運動的數(shù)值研究提供一個較為適用和有效的手段。

1 數(shù)值模型

1.1 SPH方法簡介

光滑粒子流體動力學(SPH)是一種拉格朗日無網(wǎng)格方法,使用一組粒子離散化連續(xù)介質(zhì)，用于流體動力學模擬時，SPH根據(jù)每個粒子周圍的物理特性，在每個粒子的位置對離散的Navier-Stokes方程進行局部積分。每個粒子有一個相鄰粒子的集合，基于距離的函數(shù)決定，相關(guān)的特征長度或平滑長度通常用h表示。通過在每個時間步長計算出每個粒子新的物理量，然后每個粒子根據(jù)更新的值發(fā)生移動。

SPH使用基于插值函數(shù)的積分方程將連續(xù)流體動力學的守恒定律從其偏微分形式轉(zhuǎn)化為適合于基于粒子模擬的形式。SPH函數(shù)核近似的表達式為：

(1)

式中，r為粒子的矢徑；下標a表示流體粒子；下標b為該粒子的相鄰粒子；mb為相鄰粒子的質(zhì)量;ρb為相鄰粒子的密度;W(r,h)為核函數(shù)，表達式為：

(2)

1.2 控制方程

粒子移動滿足連續(xù)介質(zhì)的動量守恒方程為：

(3)

式中，v為粒子速度矢量；ρ為粒子密度;p為壓力梯度;g為重力加速度;Γ為粘性耗散項;這里采用層流粘度和亞粒子尺度(Sub-Particle Scale，SPS)湍流模型[14]。

SPH求解流體運動時整個流場的密度是弱可壓縮的，即每個粒子的質(zhì)量保持不變，但粒子相關(guān)的密度會產(chǎn)生波動。這些密度變化是通過求解連續(xù)性方程來計算的：

(4)

式中，vab為粒子間的速度矢量差。

在流體被視為弱可壓縮性時，需要引入狀態(tài)方程來確定基于粒子密度的流體壓力。壓力和密度之間的關(guān)系遵循以下表達式：

(5)

時間步長控制方程采用可變步長格式，滿足(Courant-Friedrichs-Lewy,CFL)條件，表達式為：

(6)

式中，|fa|為單位質(zhì)量力;cs為瞬時聲速;η為阻尼系數(shù);CCFL=0.2為CFL常數(shù)。

浮體運動是通過浮體與流體粒子的相互作用力來驅(qū)動，相互作用力的求解主要通過以下步驟實現(xiàn)： ①假設浮體是剛性的；②根據(jù)指定的核函數(shù)和平滑長度，判定浮體邊界粒子所有的相鄰流體粒子；③計算作用在每個邊界粒子上的凈力并累積之和為浮體所受的作用力。每個邊界粒子k承受的單位質(zhì)量力為：

fk=∑fka,

(7)

式中，fka為流體粒子a在邊界粒子k上施加的單位質(zhì)量的力。

使用剛體動力學的基本方程式求解浮體的運動，表達式為：

(8)

(9)

式中，M為浮體的質(zhì)量；V為浮體的速度矢量；mk為邊界粒子k的質(zhì)量;I為浮體的轉(zhuǎn)動慣量;ω為浮體轉(zhuǎn)動的角速度矢;rk為邊界粒子k的矢徑;R0為浮體的質(zhì)心。

在求得每一時間步的浮體運動的速度和角速度矢基礎上，還需對浮體的邊界粒子k的速度矢量進行更新，計算式為：

vk=V+ω×(rk-R0),

(10)

2 算例驗證

2.1 活塞式造波池中的波浪傳播

圖1 活塞式造波池計算模型

造波板往復運動計算式為：

(11)

式中，H為波高；T為周期；d為靜水深度；L為波長;S0為活塞行程;e(t)為造波板t時刻的位移。

生成的二階斯托克斯波水位計算式為：

(12)

式中，x為水平方向坐標；η(x,t)為t時刻的水位。

波浪水位SPH數(shù)值計算結(jié)果與二階斯托克斯波解析解的對比如圖2所示。圖2給出了活塞行程S0=0.165 832 m，波高H=0.15 m，周期T=2.0 s，波長L=4.523 67 m條件下，x=6.0 m處波浪水位數(shù)值解與斯托克斯波解析解的對比結(jié)果。從圖2中可以看出，波浪水位的數(shù)值解與斯托克斯波水位解析解在4 s之后重合度較高。進一步提取速度測點(x=6.0 m，z=0.4 m)x與z方向速度并與斯托克斯波速度解析解進行對比，結(jié)果如圖3所示。由圖3可知，測點x與z方向速度數(shù)值解與斯托克斯波速度解析解4 s之后同樣高度吻合。需要說明的是在圖2和圖3中，前4 s數(shù)值解與解析解存在明顯差異，原因是數(shù)值模擬時，波浪從生成到傳播至測點有一段距離，在這段時間內(nèi)測點水位和速度變化幅度很小，等波浪傳播過來才會產(chǎn)生明顯水位波動，進而產(chǎn)生相應的x與z方向速度。綜合以上說明，研究數(shù)值模型數(shù)值模擬活塞式造波池中的波浪傳播過程精度較高。

圖2 波浪水位SPH數(shù)值計算結(jié)果與二階斯托克斯波解析解的對比

圖3 測點x與z方向速度SPH數(shù)值計算結(jié)果與二階斯托克斯波解析解的對比

2.2 水平圓柱自由跌落入水

以水平圓柱自由跌落入水作為算例來驗證數(shù)值模型模擬流體驅(qū)動浮體運動的正確性和適用性。水平圓柱自由跌落入水的計算模型如圖4所示。計算域x方向長10.0 m，z方向高16.0 m，靜水深度為14.0 m，水的密度ρl=1 000 kg/m3。水平圓柱直徑D=1.0 m，密度ρs=1 200 kg/m3。z方向重力加速度大小為9.81m/s2。初始時刻圓柱質(zhì)心與靜水面高差為0，釋放后圓柱穿過自由液面完成入水。

圖4 水平圓柱自由跌落入水計算模型示意圖

以自由液面為z方向坐標原點，分別提取水平圓柱質(zhì)心z方向位移和速度隨時間的變化曲線，如圖5和圖6所示。從圖5、圖6中可以看出，由于圓柱密度大于水的密度，圓柱入水過程中不斷下降，質(zhì)心z方向位移隨時間增加而一直呈現(xiàn)減小趨勢。反觀圓柱質(zhì)心z方向速度隨時間的變化曲線，呈現(xiàn)先減小后些許增加再些許下降的振蕩趨勢。將水平圓柱質(zhì)心z方向位移和速度的SPH數(shù)值模型計算結(jié)果與Fekken[15]的數(shù)值計算結(jié)果對比可以看出，位移與速度曲線的吻合度均較高，說明本研究中數(shù)值模型對于流體驅(qū)動浮體運動的模擬具有較高的精度和適用性。

圖5 水平圓柱質(zhì)心z方向位移隨時間的變化曲線 圖6 水平圓柱質(zhì)心z方向速度隨時間的變化曲線

3 活塞式造波池中箱式浮體的運動

活塞式造波池中箱式浮體運動計算模型如圖7所示。計算模型的幾何尺寸為長×寬×高=16.0 m×6.0 m×8.0 m?；钊皆觳ò宄叽鐬殚L×寬×高=0.3 m×6.0 m×7.0 m，離計算模型左邊界0.7 m。造波板行程為1.0 m，其水平往復周期運動過程為先向x正方向運動0.5 m后返回原位置再向x負方向運動0.5 m再回到原位置，其運動周期為10/3 s。正方體箱式浮體位于離左邊界10.0 m的水池中央，尺寸為長×寬×高=2.0 m×2.0 m×2.0 m，密度為ρs=1 250 kg/m3。造波池中的靜水深度為4.0 m，水的密度為ρl=1 000 kg/m3。z方向重力加速度大小為9.81 m/s2。初始時刻計算模型的流場速度為0，箱式浮體運動時長為16.7 s。

圖7 活塞式造波池中箱式浮體運動計算模型及幾何尺寸

為對波浪驅(qū)動下箱式浮體運動做深入剖析，現(xiàn)將箱式浮體六自由度運動描繪出來，如圖8所示。其中，沿3個正交坐標軸的平移運動包括垂向、橫向和縱向運動；繞3個正交坐標軸的旋轉(zhuǎn)運動包括偏航、俯仰和橫搖轉(zhuǎn)動。進一步，將箱式浮體六自由度運動隨時間變化的數(shù)據(jù)進行，位移和轉(zhuǎn)動結(jié)果分別如圖9和圖10所示。從圖9中可以看出，在波浪驅(qū)動下，箱式浮體縱向和垂向位移曲線均呈現(xiàn)規(guī)律的波動狀態(tài)，說明浮體位移對波浪的響應以縱向和垂向為主，橫向運動受波浪驅(qū)動的影響較小。由于計算模型右邊界沒有進行消波處理，活塞持續(xù)水平往復運動產(chǎn)生的波浪激勵逐步增強，使得箱式浮體縱向位移振幅呈現(xiàn)明顯增加趨勢。而垂向的浮力作用主要由密度差決定，所以箱式浮體垂向位移振幅的增加趨勢并不明顯。箱式浮體繞正交坐標軸的旋轉(zhuǎn)運動角度曲線如圖10所示。由圖10中可以看出，繞y軸俯仰轉(zhuǎn)動以及繞x軸橫搖轉(zhuǎn)動的角度曲線均呈現(xiàn)規(guī)律的波動狀態(tài)，說明箱式浮體轉(zhuǎn)動對波浪的響應以俯仰轉(zhuǎn)動和橫搖轉(zhuǎn)動為主，偏航轉(zhuǎn)動影響受波浪驅(qū)動的影響較小。同樣，受波浪激勵逐步增強的影響，箱式浮體繞y軸俯仰轉(zhuǎn)動的角度振幅呈現(xiàn)逐漸增強趨勢，繞x軸橫搖轉(zhuǎn)動角度振幅的增加趨勢并不明顯。

圖8 箱式浮體六自由度運動示意圖

圖9 箱式浮體沿正交坐標軸平移運動位移 圖10 箱式浮體繞正交坐標軸旋轉(zhuǎn)運動角度

4 結(jié)論

研究通過SPH方法建立了一個用于研究活塞式造波池中浮體運動的數(shù)值模型，并通過活塞式造波池中的波浪傳播以及水平圓柱自由跌落入水分別驗證了該數(shù)值模型在活塞式造波以及流體驅(qū)動浮體運動兩個方面的正確性和適用性。在此基礎上，對活塞式造波池中箱式浮體的運動進行了數(shù)值研究，并對箱式浮體六自由度運動的數(shù)據(jù)和響應規(guī)律進行了深入剖析，得到了以下結(jié)論：箱式浮體位移對波浪的響應以縱向和垂向為主，橫向運動受波浪驅(qū)動的影響較??；箱式浮體縱向位移振幅易受波浪激勵的影響，而垂向位移振幅受波浪激勵的影響較小。箱式浮體轉(zhuǎn)動對波浪的響應以俯仰轉(zhuǎn)動和橫搖轉(zhuǎn)動為主，偏航轉(zhuǎn)動受波浪驅(qū)動的影響較??；箱式浮體繞y軸俯仰轉(zhuǎn)動的角度振幅易受波浪激勵的影響，而繞x軸橫搖轉(zhuǎn)動角度振幅受波浪激勵的影響較小。