李增亮, 王樂峰, 董祥偉, 杜明超, 荊正軍, 王雨婷

(1.中國石油大學(xué)(華東)機電工程學(xué)院,山東青島 266580; 2.中國石油大學(xué)(華東)石油工業(yè)訓(xùn)練中心,山東青島 266580;3.中國石油海洋工程有限公司青島分公司，山東青島 266555)

波浪與浮體的相互作用現(xiàn)象廣泛存在于海洋能源領(lǐng)域,例如海洋平臺在波浪中的運動響應(yīng)[1],浮子式波浪能發(fā)電裝置的波浪能轉(zhuǎn)換過程[2]等。近年來,國內(nèi)外在波浪能轉(zhuǎn)換方面做了很多研究,世界上第一個關(guān)于波浪能發(fā)電技術(shù)的專利誕生于1799年[3],日本于20世紀80年代初建成了總裝機高達1 250 kW的波浪能轉(zhuǎn)換裝置[4],中國的廣州能源研究所于1984年研制成功的波浪能轉(zhuǎn)換裝置已在沿海海域大規(guī)模投入使用[5]。波浪與浮體的相互作用是典型的流固耦合問題,同時涉及自由表面流動、浮體結(jié)構(gòu)物的運動等。學(xué)者們采用數(shù)值模擬手段對波浪的水動力學(xué)特性進行了研究。王永學(xué)[6]利用VOF方法建立無反射造波數(shù)值波浪水槽；董志等[7]利用VOF方法建立數(shù)值波浪水槽并對造波消波方法進行研究；胡杭輝等[8]進行了基于Fluent的二維非線性數(shù)值波浪水槽構(gòu)造及驗證的研究；臧志鵬等[9]利用試驗裝置研究波浪發(fā)電系統(tǒng)振蕩浮體的運動特性。但是,當(dāng)涉及波浪與浮體相互作用的數(shù)值模擬時,由于網(wǎng)格特性的限制和浮體的拉格朗日運動特性,傳統(tǒng)的基于歐拉方法(如有限體積法)的數(shù)值模擬手段并不適用于該類問題的求解。近年來,無網(wǎng)格粒子法逐漸興起,光滑粒子動力學(xué)方法(SPH)是最具代表性的方法之一[10]。倪興也等[11]采用SPH方法建立二維數(shù)值波浪水槽,模擬了推板造波過程;常江[12]采用SPH方法建立三維數(shù)值波浪水槽,得到的波浪運動特性與試驗結(jié)果吻合較好。作為一種粒子法,SPH方法的主要優(yōu)點在于不需要畫網(wǎng)格,在處理自由表面流動與固體的流固耦合過程中較有優(yōu)勢。筆者將SPH方法應(yīng)用于波浪與浮體相互作用的數(shù)值仿真研究,借助開源程序DualSPHysics建立三維數(shù)值波浪水槽,并搭建造波水槽試驗臺,通過試驗對比驗證數(shù)值模型的適用性。

1 數(shù)值方法及模型

1.1 SPH方法原理

光滑粒子動力學(xué)方法(smoothed particle hydrodynamics,SPH)是一種拉格朗日無網(wǎng)格法。與有限體積等網(wǎng)格法不同,SPH方法利用一組粒子離散連續(xù)體計算域。在SPH流體力學(xué)仿真中,根據(jù)周圍顆粒物的物理性質(zhì),Navier Stokes方程在每個粒子處進行局部積分,積分域內(nèi)包含粒子的集合是由基于距離的函數(shù)來確定的。例如,對于二維問題,該積分域是以相關(guān)粒子為圓心的圓形區(qū)域;對于三維問題,積分域是以相關(guān)粒子為球心的球形區(qū)域。在SPH方法中,物理量賦值在每一個粒子上,在每一個時間步粒子各物理量都會被重新計算,然后粒子根據(jù)這些重新計算的數(shù)據(jù)產(chǎn)生新的位移。因此,SPH方法是一種拉格朗日數(shù)值模擬方法。

假定一任意場函數(shù)F(r),其可以通過核函數(shù)近似寫成積分表達形式:

(1)

式中,W(r-r′,h)為核函數(shù);h為核函數(shù)的光滑長度。

式(1)又稱為場函數(shù)F的核近似表達。由于計算域被離散成一系列粒子,F函數(shù)的粒子近似可寫成:

(2)

式中,下標(biāo)a代表當(dāng)前SPH粒子;下標(biāo)b代表粒子a的鄰域粒子(b粒子位于a粒子的支持域內(nèi));mb為粒子b的質(zhì)量;ρb為粒子b的密度。

1.2 核函數(shù)選擇

在SPH方法中,存在多種形式的核函數(shù),核函數(shù)的選取將直接影響數(shù)值計算的精度。DualSPHysics程序中內(nèi)置了2種核函數(shù),分別為五次樣條核函數(shù)和三次樣條核函數(shù)。經(jīng)過對比發(fā)現(xiàn),在其他條件不變的情況下,五次樣條核函數(shù)得到的波浪流場的壓力分布結(jié)果要優(yōu)于三次樣條核函數(shù),因此選取了五次樣條核函數(shù)[13]:

(3)

其中

式中,q為無量綱距離;r為任意兩粒子a與b之間的實際距離;h為光滑長度。

對于任意粒子a,其領(lǐng)域粒子是由粒子a的光滑半徑?jīng)Q定的,光滑半徑又與光滑長度h有關(guān)。對于式(3),當(dāng)q大于2時核函數(shù)的值為0;只有q小于2時,核函數(shù)才是有效的,這也說明SPH核函數(shù)是緊致的。

1.3 流體控制方程

流體運動由流體控制方程描述,即納維斯托克斯方程(N-S),包括質(zhì)量守恒方程和動量守恒方程。質(zhì)量守恒方程為

(4)

式中,ρ為流體密度,kg/m3;v為流體速度,m/s。

在SPH中,方程(4)可以寫成粒子離散形式,表示為

(5)

其中

vab=va-vb.

式中,mb為粒子b質(zhì)量,kg;vab為粒子速度,m/s;Wab為a和b粒子之間核函數(shù)值。

圖1 光滑核函數(shù)及其導(dǎo)數(shù)圖像Fig.1 Smooth kernel functions and derivatives of smooth kernel functions

動量守恒方程為

(6)

式中,t為時間,s;p為壓力,Pa;g為自由落體加速度,m/s2;Γ為耗散項。

采用粒子近似方法,方程(6)可以寫成SPH離散形式:

(7)

式中,pa和pb分別為粒子a、b處壓力;ρa和ρb分別為粒子a、b處密度。

Πab為人工黏性力[14],其作用是消除數(shù)值計算時計算域中的非物理震蕩。人工黏性力為

(8)

其中

式中,v為粒子速度,m/s;α為常數(shù),在自由表面流動模擬中一般取0.01[15];Cab為聲速的平均值,即a粒子和b粒子代表的流體的聲速平均(流體為水介質(zhì),因此a粒子和b粒子的聲速相同)。

在本文模擬中,流體被當(dāng)作弱可壓縮的,采用狀態(tài)方程來計算流體壓力。狀態(tài)方程表達為流體壓力和密度之間的關(guān)系式[16]:

(9)

其中

式中,ρ0為參考密度,取為1 000 kg/m3;c0為流體聲速。

水的物理聲速為1 500 m/s,但在實際模擬過程中,由于流體是弱可壓縮性的,一般可通過給定聲速一個較低值,只要保證密度的波動低于1%即可;實際的聲速可計算為c=10vmax,其中vmax為計算域中流體的最大速度。

1.4 浮體受力及運動方程

浮體與流體相互作用過程中,將浮體考慮為剛體,并也離散成SPH粒子。根據(jù)積分步驟,可以計算出每一邊界粒子k受到的周圍流體粒子所給出的力fk:

(10)

mkfka=-mafak.

(11)

式中,fka為流體粒子a作用在邊界粒子k上的作用力(a粒子位于k粒子的支持域之內(nèi)),N;fak為邊界粒子k作用在流體粒子a上的反作用力,N;ma和mk分別為流體粒子和邊界粒子的質(zhì)量,kg。

圖2為浮體的邊界條件。

圖2 浮體的邊界條件Fig.2 Boundary conditions for floating bodies

對于運動的剛體,可通過剛體的運動方程求出,包括平動方程和轉(zhuǎn)動方程:

(12)

(13)

式中,M為剛體質(zhì)量,kg;I為慣性矩,m4;Ω為角速度,rad/s;R0代表質(zhì)心。

對式(12)、(13)進行積分即可得到剛體在某時間的速度和角速度。由于剛體的SPH粒子是附著在剛體上的,可以認為是質(zhì)點,則每一邊界粒子的速度uk為

uk=v+Ω(rk-R0).

(14)

1.5 數(shù)值模型

數(shù)值模擬是基于DualSPHysics完成的。DualSPHysics是一款開源的SPH程序,由c++語言編寫,借助于FreeCAD前處理接口程序可以實現(xiàn)在Windows平臺對流體力學(xué)問題進行數(shù)值建模和SPH模擬。依據(jù)實驗室內(nèi)波浪水槽的實際尺寸,建立波浪和浮體相互作用的數(shù)值模型,如圖3所示。模型包含造波板、浮體、水、消波板及水槽。其中造波板、浮體、消波板及水槽壁面被處理為剛性邊界。給定造波板一定的運動規(guī)律,推擠流體產(chǎn)生波浪運動行為。浮體被建模為可以自由運動的剛體,浮體與水的流固耦合作用由動態(tài)邊界法(dynamic boundary)實現(xiàn)。圖4為計算域示意圖。計算域離散為一系列初始均勻分布的SPH粒子,粒子間距取值為8 mm。

圖3 三維數(shù)值模型Fig.3 Three-dimensional numerical model

圖4 離散后的三維數(shù)值模型Fig.4 Discrete three-dimensional numerical model

DualSPHysics中內(nèi)置了多種時間積分格式,選用辛格式(symplectic scheme)。在第一階段,其位移和密度可表示為

(15)

(16)

(17)

其中

(18)

(19)

(20)

仿真結(jié)果以數(shù)據(jù)表的形式輸出,借助于第三方后處理軟件Paraview完成結(jié)果的可視化。仿真過程中設(shè)置輸出數(shù)據(jù)頻率為每0.01 s輸出一次,仿真的總物理時間為25 s。

造波采用推板式造波,水槽右側(cè)采用斜坡式消波。水槽整體長度為3 m,高度為45 cm,寬度為50 cm;造波板寬度為492 mm,高度為400 mm,厚度為1.4 mm,造波板與水槽底部和水槽兩側(cè)的距離均為4 mm;浮體距離水槽左端80 cm且距離水槽兩個側(cè)面的距離相等,使得浮體處于有效工作區(qū)域。

為了研究浮體形狀的影響,選取扁平型、球型和圓柱型3種浮體。首先研究圓柱型浮體,圓柱半徑為50 mm,高度分別為30、50和70 mm,如圖5所示;球型浮體模型,半徑為45.428 mm;扁平型的浮體由兩個球缺疊加在一起,高度為50 mm,球缺半徑為108.33 mm;3種浮體如圖6所示。浮體材料的密度設(shè)置為600 kg/m3,仿真過程中粒子間距為8 mm,水深100 mm。

圖5 不同高度圓柱型浮體Fig.5 Cylindrical floater with different heights

圖6 不同形狀浮體外觀Fig.6 Floater of different shapes

1.6 數(shù)值模型驗證

根據(jù)已建立好的數(shù)值模型進行造波模擬。波形采用二階斯托克斯波,表述為

(21)

式中,S0為造波板行程,m;ω為角速度,rad/s;δ為初相,rad。

(22)

其中

k=2π/L,

式中,H為波高,m;d為水深,m;L為波長,m。

二階斯托克斯波為

e(t)=e1(t)+e2(t).

(23)

設(shè)置波高0.03 m、周期1 s、水深0.1 m。圖7為不同時刻的波形圖,圖8為采用SPH方法模擬得到的波形與理論解對比。圖8中SPH模擬得到的波形與理論解吻合較好,說明所建立的SPH模型能夠有效預(yù)測推板造波水槽中波浪的自由表面流動特性。

圖7 不同時刻波形Fig.7 Different instants of simulation with regular waves

圖8 SPH和理論解波形對比Fig.8 Comparison of SPH and theoretical wave surface elevation

2 波浪水槽試驗

2.1 試驗裝置和方案

搭建造波水槽試驗臺如圖9所示。利用該試驗裝置進行波浪與浮體相互作用試驗,將試驗結(jié)果與仿真結(jié)果進行對比。試驗臺由造波機構(gòu)和水槽兩部分組成。造波機構(gòu)包括電源、伺服電機控制器、伺服電機驅(qū)動器、750 W伺服電機、聯(lián)軸器、絲杠導(dǎo)軌、造波板,如圖10所示。水槽包括浮體、消波板、玻璃容器和支架,如圖11所示。造波機構(gòu)的工作原理為:伺服電機控制器控制電機旋轉(zhuǎn),絲杠導(dǎo)軌將電機的旋轉(zhuǎn)運動轉(zhuǎn)變?yōu)樵觳ò宓耐鶑?fù)直線運動,造波板推擠流體產(chǎn)生波浪。水槽用透明玻璃制成,用高速攝像機記錄浮體在波浪中的運動過程。后續(xù)通過圖像處理,記錄浮體的位置,通過描點的方式繪出浮體的運動軌跡,從而將試驗數(shù)據(jù)與仿真數(shù)據(jù)進行對比。

圖9 波浪水槽試驗臺Fig.9 Wave-making experimental device

圖10 造波機構(gòu)示意圖Fig.10 Sketch map of wave making mechanism

圖11 水槽部分示意圖Fig.11 Sketch map of flume part

造波板的運動速度與電機的轉(zhuǎn)速關(guān)系為

v′=kn.

(24)

式中,v′為造波板的運動速度,mm/min;n為電機轉(zhuǎn)速,r/min;k為絲杠導(dǎo)程,本裝置中k為4 mm。

2.2 水槽中浮體運動試驗及仿真

圖12為浮體初始時刻在流體中保持靜止時的位置。根據(jù)阿基米德浮力定律可算出浮體在液體中靜止時的吃水深度。

Ff=ρlgVd.

(25)

式中,Ff為浮體受到的浮力,N;ρl為液體密度,此處為水的密度,1 000 kg/m3;Vd為浮體排開液體的體積,m3。

圖12 浮體在液面中靜止時試驗與仿真對比Fig.12 Comparison of experiment and simulation of floater at rest in liquid surface

圖13為不同時刻試驗與仿真對比。造波板運動規(guī)律如表1所示。為了保證仿真與試驗條件一致,仿真過程中造波板運動規(guī)律與試驗相同。

圖13 仿真與試驗不同時刻對比Fig.13 Comparison of simulation and experiment at different time

表1 試驗與仿真過程造波板運動規(guī)律(部分數(shù)據(jù))
Table1Motionlawofwave-makingplateinexperimentsandsimulations(some data)

序號速度/(m·s-1)時間/s行程/m 10.16660.360.06 2-0.16660.360.06 30.16660.360.06 4-0.16660.360.06 ………… 350.16660.360.06 36-0.16660.360.06

圖14為水平方向和豎直方向浮體位置隨時間的變化規(guī)律。由圖14可知,若忽略裝置造波板剛性不足的影響以及試驗后處理過程中對浮體位置進行測量時誤差的影響,試驗與仿真結(jié)果基本一致,可認為采用本文中的仿真方法能夠有效預(yù)測浮體在波浪環(huán)境下的運動規(guī)律,進而可以通過仿真來研究其他因素對波浪與浮體相互作用產(chǎn)生的影響。

圖14 浮體水平及豎直方向隨時間運動規(guī)律試驗與仿真對比Fig.14 Comparison of experiment and simulation on law of horizontal and vertical motion of floating body with time

3 浮體在波浪水槽中的運動

3.1 浮體形狀影響

仿真過程中為了保證各個浮體只有形狀的區(qū)別,浮體的密度、體積均相同,波浪參數(shù)相同。

仿真過程中造波板的運動規(guī)律如表2所示,圖15為不同形狀的浮體在水平和豎直方向浮體位置隨時間運動規(guī)律。

表2 仿真過程造波板運動規(guī)律(部分數(shù)據(jù))

通過觀察圖15可以看出,水平方向,浮體沿著波浪方向不斷前進,初始位置相同,初期位移曲線基本吻合,后期不同形狀浮體的位移之差逐漸加大,扁平型在水平方向產(chǎn)生的位移最大,圓柱型次之,球型最小;豎直方向,初始位置相同,不同形狀的浮體振幅基本相同,后期扁平型浮體相對于其他兩種浮體縱坐標(biāo)較高,圓柱型和球型兩種浮體曲線基本吻合。說明浮體的形狀對波浪與浮體相互作用影響較大,特別是在水平方向?qū)Ω◇w的運動規(guī)律有較大影響。

圖15 不同形狀浮體水平及豎直方向隨時間運動規(guī)律Fig.15 Comparison of horizontal and vertical motion of floater with different shapes with time

圖16為浮體橫搖和縱搖隨時間變化規(guī)律。其方向如圖17所示。橫搖與縱搖的大小體現(xiàn)出浮體搖晃的程度,在浮子式波浪能發(fā)電中,應(yīng)盡量減小橫搖與縱搖,降低橫搖與縱搖對浮子式波浪能發(fā)電裝置的影響。由圖16可以看出,球型的搖晃角度最大,扁平型和圓柱型的接近。從圖13可以直觀地看出,隨著波峰與波谷交替變換,浮體必然會產(chǎn)生縱搖。由圖16可以看出,由于波浪傳播方向的影響,扁平型和圓柱型縱搖的角度要遠大于橫搖的。

圖16 不同形狀浮體橫搖及縱搖隨時間變化規(guī)律Fig.16 Variation of rolling and pitching of floater with different shapes with time

圖17 橫搖與縱搖所表示的方向Fig.17 Rolling and pitching that direction

3.2 浮體尺寸影響

圖18表示不同高度圓柱體浮體在16 s時刻水平方向發(fā)生的位移和水平方向受到的力。由于浮體高度逐漸增加,浮體質(zhì)量會逐漸增大,浮體所受到的阻力會逐漸增大。由圖18可以看出,隨著浮體高度增加,浮體在水平方向位移逐漸減小,受到的力逐漸增大。圖19為浮體隨著高度增加發(fā)生傾斜現(xiàn)象。

圖18 不同高度圓柱體浮體的位移與受力Fig.18 Displacement and force of cylindrical floater with different heights

圖19 不同高度圓柱型浮體產(chǎn)生傾斜現(xiàn)象Fig.19 Inclination of cylindrical floater at different heights

由圖19可知,高度越高,傾斜角度越大。浮子式波浪能發(fā)電裝置中,浮體傾斜產(chǎn)生的扭矩對于發(fā)電裝置非常不利,因此在發(fā)電裝置中應(yīng)該盡可能避免浮體發(fā)生傾斜現(xiàn)象,在設(shè)計浮體參數(shù)時應(yīng)使浮體在波浪中避免發(fā)生傾斜現(xiàn)象。

3.3 浮體密度影響

圖20為不同密度浮體在一個波峰經(jīng)過浮體時的越浪現(xiàn)象對比。由圖20可知,當(dāng)波峰逐漸向浮體靠近時,浮體左側(cè)被迫升高,右側(cè)降低;當(dāng)波峰達到浮體時,由于受到浮體的阻礙,使波浪質(zhì)點上升,對于密度較大的浮體,其不能立即跟隨波峰升高,從而導(dǎo)致波浪越過浮體,發(fā)生越浪現(xiàn)象。密度為600 kg/m3的浮體幾乎不會產(chǎn)生越浪現(xiàn)象,而隨著密度增加越浪現(xiàn)象越來越明顯。密度為800 kg/m3的浮體即能觀察到越浪現(xiàn)象,密度為900 kg/m3的浮體越浪效果更加明顯。圖21為有無越浪現(xiàn)象對比,圖中左側(cè)為密度為600 kg/m3的浮體,無越浪現(xiàn)象,右側(cè)為密度為900 kg/m3的浮體,有明顯的越浪現(xiàn)象。

圖20 不同密度浮體越浪效果對比Fig.20 Comparison of overtopping effects of floater with different densities

圖21 有無越浪現(xiàn)象對比Fig.21 Comparison of overtopping phenomenon and non-overtopping phenomenon

4 結(jié) 論

(1)相比于歐拉網(wǎng)格法,無網(wǎng)格SPH方法的建模流程更便捷,采用粒子直接填充流體計算域,流固耦合過程不受限于網(wǎng)格拓撲關(guān)系,更適應(yīng)復(fù)雜幾何形狀的計算域,適用于浮體的大位移運動模擬。

(2)通過波浪理論驗證了SPH模型對波浪運動波形模擬的正確性。在此基礎(chǔ)上進一步模擬了浮體與波浪作用過程,得到的浮體位移-時間關(guān)系與試驗數(shù)據(jù)吻合較好。同時,捕捉到了越浪、波浪作用下的浮體側(cè)傾等現(xiàn)象,展現(xiàn)了無網(wǎng)格SPH方法在自由表面流動和浮體相互作用模擬中的優(yōu)勢。

(3)不同形狀的浮體在波浪中的運動規(guī)律具有明顯差別,相同形狀的浮體不同的尺寸參數(shù)、密度對浮體的運動也會產(chǎn)生較大影響。