基于無網(wǎng)格法的浮體與波浪相互作用模擬與試驗(yàn)驗(yàn)證

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

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

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

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

1.1 SPH方法原理

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

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

(1)

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

式(1)又稱為場(chǎng)函數(shù)F的核近似表達(dá)。由于計(jì)算域被離散成一系列粒子,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ù)值計(jì)算的精度。DualSPHysics程序中內(nèi)置了2種核函數(shù),分別為五次樣條核函數(shù)和三次樣條核函數(shù)。經(jīng)過對(duì)比發(fā)現(xiàn),在其他條件不變的情況下,五次樣條核函數(shù)得到的波浪流場(chǎng)的壓力分布結(jié)果要優(yōu)于三次樣條核函數(shù),因此選取了五次樣條核函數(shù)[13]:

(3)

其中

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

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

1.3 流體控制方程

流體運(yùn)動(dòng)由流體控制方程描述,即納維斯托克斯方程(N-S),包括質(zhì)量守恒方程和動(dòng)量守恒方程。質(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

動(dòng)量守恒方程為

(6)

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

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

(7)

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

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

(8)

其中

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

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

(9)

其中

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

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

1.4 浮體受力及運(yùn)動(dòng)方程

浮體與流體相互作用過程中,將浮體考慮為剛體,并也離散成SPH粒子。根據(jù)積分步驟,可以計(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

對(duì)于運(yùn)動(dòng)的剛體,可通過剛體的運(yùn)動(dòng)方程求出,包括平動(dòng)方程和轉(zhuǎn)動(dòng)方程:

(12)

(13)

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

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

uk=v+Ω(rk-R0).

(14)

1.5 數(shù)值模型

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

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

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

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

(15)

(16)

(17)

其中

(18)

(19)

(20)

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

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

為了研究浮體形狀的影響,選取扁平型、球型和圓柱型3種浮體。首先研究圓柱型浮體,圓柱半徑為50 mm,高度分別為30、50和70 mm,如圖5所示;球型浮體模型,半徑為45.428 mm;扁平型的浮體由兩個(gè)球缺疊加在一起,高度為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ù)值模型驗(yàn)證

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

(21)

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

(22)

其中

k=2π/L,

式中,H為波高,m;d為水深,m;L為波長(zhǎng),m。

二階斯托克斯波為

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

(23)

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

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

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

2 波浪水槽試驗(yàn)

2.1 試驗(yàn)裝置和方案

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

圖9 波浪水槽試驗(yàn)臺(tái)Fig.9 Wave-making experimental device

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

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

造波板的運(yùn)動(dòng)速度與電機(jī)的轉(zhuǎn)速關(guān)系為

v′=kn.

(24)

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

2.2 水槽中浮體運(yùn)動(dòng)試驗(yàn)及仿真

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

Ff=ρlgVd.

(25)

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

圖12 浮體在液面中靜止時(shí)試驗(yàn)與仿真對(duì)比Fig.12 Comparison of experiment and simulation of floater at rest in liquid surface

圖13為不同時(shí)刻試驗(yàn)與仿真對(duì)比。造波板運(yùn)動(dòng)規(guī)律如表1所示。為了保證仿真與試驗(yàn)條件一致,仿真過程中造波板運(yùn)動(dòng)規(guī)律與試驗(yàn)相同。

圖13 仿真與試驗(yàn)不同時(shí)刻對(duì)比Fig.13 Comparison of simulation and experiment at different time

表1 試驗(yàn)與仿真過程造波板運(yùn)動(dòng)規(guī)律(部分?jǐn)?shù)據(jù))
Table1Motionlawofwave-makingplateinexperimentsandsimulations(some data)

序號(hào)速度/(m·s-1)時(shí)間/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為水平方向和豎直方向浮體位置隨時(shí)間的變化規(guī)律。由圖14可知,若忽略裝置造波板剛性不足的影響以及試驗(yàn)后處理過程中對(duì)浮體位置進(jìn)行測(cè)量時(shí)誤差的影響,試驗(yàn)與仿真結(jié)果基本一致,可認(rèn)為采用本文中的仿真方法能夠有效預(yù)測(cè)浮體在波浪環(huán)境下的運(yùn)動(dòng)規(guī)律,進(jìn)而可以通過仿真來研究其他因素對(duì)波浪與浮體相互作用產(chǎn)生的影響。

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

3 浮體在波浪水槽中的運(yùn)動(dòng)

3.1 浮體形狀影響

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

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

表2 仿真過程造波板運(yùn)動(dòng)規(guī)律(部分?jǐn)?shù)據(jù))

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

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

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

圖16 不同形狀浮體橫搖及縱搖隨時(shí)間變化規(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時(shí)刻水平方向發(fā)生的位移和水平方向受到的力。由于浮體高度逐漸增加,浮體質(zhì)量會(huì)逐漸增大,浮體所受到的阻力會(huì)逐漸增大。由圖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)生的扭矩對(duì)于發(fā)電裝置非常不利,因此在發(fā)電裝置中應(yīng)該盡可能避免浮體發(fā)生傾斜現(xiàn)象,在設(shè)計(jì)浮體參數(shù)時(shí)應(yīng)使浮體在波浪中避免發(fā)生傾斜現(xiàn)象。

3.3 浮體密度影響

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

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

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

4 結(jié) 論

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

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

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