振蕩浮子波能裝置波浪沖擊特性數(shù)值模擬研究

奚泉，唐漁瀅，于通順，于春明，曾興井，徐琨

1.中國電建集團昆明勘測設計研究院有限公司，云南 昆明 650051

2.青島市交通工程質(zhì)量安全監(jiān)督站，山東 青島 266101

3.中國海洋大學 工程學院，山東 青島 266400

4.中國葛洲壩集團市政工程有限公司，湖北 宜昌 443000

在“雙碳”的戰(zhàn)略目標下，海洋可再生能源迎來發(fā)展新格局[1]。波浪能因具有開發(fā)難度較低、獲取方式簡單等優(yōu)勢，被認為是品質(zhì)最高的海洋能[2-3]。振蕩浮子波能發(fā)電裝置[4-5]因結構簡單、能量轉(zhuǎn)換效率高等優(yōu)勢是目前常見的三大波浪能利用技術之一。

振蕩浮子裝置大多是由波能俘獲裝置、質(zhì)量體、系泊系統(tǒng)以及發(fā)電裝置組成，其工作原理是利用俘獲裝置在波浪的作用下振蕩起伏來達到汲取波浪能的目的，之后再由能量攝?。╬ower takeoff，PTO）系統(tǒng)實現(xiàn)波浪能向電能轉(zhuǎn)化。振蕩浮子式波浪能發(fā)電裝置的水動力性能[6]、能量轉(zhuǎn)換效率[7]及性能優(yōu)化[8]成為目前研究的熱點。美國Ocean Power Technologies (OPT)公司研發(fā)了PowerBuoy 振蕩浮子波能裝置，并對其進行了設計優(yōu)化[9]；Shi[10]和Han等[11]對振蕩浮子獲能特性進行了理論分析，并進一步提出了獲能譜的概念；于通順等[12-13]探究了振蕩浮子前側波浪爬升特性，為后續(xù)浮子豎向尺寸設計提供理論依據(jù)；Agamloh等[14]利用有限體積法分別模擬單個浮子和雙浮子在不同規(guī)則波作用下的運動過程，探究了浮子的能量輸出功率。吳必軍等[15]和余海濤等[16]均對單浮子波能裝置的俘獲寬度比進行了試驗研究，并指出波浪頻率是影響裝置俘獲寬度比的主要因素。Waters等[17]對驅(qū)動線性發(fā)電機的浮子進行了試驗研究，提出非線性外負載是影響裝置俘獲寬度比的主要因素。

通過以上研究發(fā)現(xiàn)，前人的研究中尚未開展浮子垂蕩過程中表面所受的沖擊作用的研究。振蕩浮子在復雜海況中通過大幅垂蕩獲能，垂蕩過程中浮子與波浪作用過程劇烈，忽略了浮子表面受到的強烈波浪沖擊作用，將會增大浮子長期運行的風險。因此本文基于FLOW 3D 建立波浪與浮子耦合作用的數(shù)值模型，開展浮子所受沖擊作用過程的研究，探究浮子所受沖擊荷載的作用規(guī)律。

1 數(shù)值計算方法

1.1 控制方程

對于波浪和結構物相互作用的問題，波浪被認為是不可壓縮黏性流體的運動?？紤]湍流流動，將連續(xù)性方程和雷諾平均納維-斯托克斯（Reynolds-averaged Navier-Stokes，RANS）方程作為流體運動的控制方程。假設右手笛卡爾坐標系O-xyz以未受擾動自由面的原點定義，則連續(xù)性方程可以寫作：

動量方程為

式中：ρ為流體密度，P為流體內(nèi)部壓強，VF為可流動的體積分數(shù)，u、v、w分別為在x、y、z這3 個方向上的流速分量，Ax、Ay、Az分別為x、y、z這3 個方向可流動的面積分數(shù)，Gx、Gy、Gz分別為對應x、y、z這3個方向的重力加速度，fx、fy、fz分別對應3個方向的黏滯力加速度。

1.2 紊流模型

本文采用重組化RNGk-ε模型，其表達式為

式中：PT為由剪切效應引起的湍流動能的生成項，KT,Diff和 εT,Diff為擴散項，CDIS,1和CDIS,2為模型系數(shù)，k為由流速梯度引起的紊動動能的產(chǎn)生項。

1.3 運動與碰撞模型

Flow-3D 中運動與碰撞模型可以模擬流固耦合下結構物的運動，該運動可以是指定的運動，也可以是與流體動態(tài)耦合的運動，移動物體都可以具有6 個自由度。該模型還允許存在多個（最多100 個）線性彈簧和扭力彈簧、彈性繩索以及系泊纜繩，它們可以附接到移動的物體上，并可對其施加力或扭矩。根據(jù)剛體運動學，剛體的一般運動可分為平移和旋轉(zhuǎn)運動。結構物上任何點的速度等于結構物上任意選擇的基點速度加上繞該基點旋轉(zhuǎn)引起的速度。選擇S為剛體上的基點，R為剛體上任意點，則R點速度表達式為

式中：VS為基點的速度；rR/S為點R相對于基點S的矢量距離；ω是剛體的角速度矢量，與基點的選擇無關。

1.4 自由表面追蹤的流體體積方法

FLOW-3D 中采用了流體體積方法進行自由追蹤。定義流體計算域內(nèi)獨個單元的流體體積函數(shù)為F（x,y,z,t)。F為流體所占單元體的比值，若單元全部充滿流體，則F值為1；若不存在流體時，F(xiàn)值則為0；而當F介于0～1 時，則表明該網(wǎng)格內(nèi)存在自由表面。F函數(shù)表示單位體積，并滿足方程：

2 數(shù)值模型的建立

2.1 三維數(shù)值波浪水槽的建立

建立數(shù)值波浪水槽的側面示意如圖1（a）所示，水槽長度取6 倍左右波長，寬度取浮子寬度的5 倍。為便于后續(xù)進行數(shù)值模型的驗證，選用Yu等[12]試驗中浮子的尺寸進行模型的建立：圓柱型浮子直徑為0.4 m，高度為0.4 m，吃水深度為0.12 m。為了合理描述浮子的形狀，該水槽由3 個網(wǎng)格塊組成，中間網(wǎng)格塊為浮子運動加密區(qū)，浮子放置在距離入口的2 倍波長(L)范圍內(nèi)。水槽的邊界條件定義為：最左端是可以實現(xiàn)函數(shù)造波功能的波浪，最右端為出流，水槽的底部設為壁面，上部為壓力，水槽順波向兩側對稱。如圖1（b）和圖1（c）所示。水槽尾端部分（約1/4 水槽長度）設有阻尼消波區(qū)用來抑制波浪反射，左端和右端的阻尼系數(shù)分別設為1 和8。

圖1 三維數(shù)值波浪水槽邊界設置及網(wǎng)格塊劃分

2.2 網(wǎng)格的劃分及網(wǎng)格獨立性驗證

本模型采用波高范圍內(nèi)3 種不同的網(wǎng)格數(shù)量進行了波面網(wǎng)格獨立性試驗，即設置了網(wǎng)格數(shù)量分別為10、20 和30 的3 個不同算例，流域內(nèi)網(wǎng)格的大小和數(shù)量如表1 所示。

表1 無浮子時計算流域網(wǎng)格數(shù)量和尺寸

圖2 比較了3 種網(wǎng)格尺寸下計算區(qū)域內(nèi)相同位置（2L處）自由面高程的歷時曲線，結果發(fā)現(xiàn)，加密網(wǎng)格數(shù)為20 與30 的計算結果差距較小，加密網(wǎng)格數(shù)為10 的計算結果與其他2 種網(wǎng)格數(shù)量的計算結果差距較大，但是加密網(wǎng)格數(shù)為30 的運行時間明顯增加。綜合考慮模型計算時間和精度，本文在Z方向的波高(H)范圍內(nèi)網(wǎng)格尺寸建議設置為 Δz≤H/20。

圖2 不同網(wǎng)格尺寸的自由面高程歷時曲線

為了合理描述浮子的形狀，在浮子運動加密區(qū)采用3 種不同總網(wǎng)格數(shù)量，浮子周圍網(wǎng)格數(shù)量如 表2所示，表中Δx、Δy和Δz分別為沿x、y和z方向的網(wǎng)格尺寸。所劃分網(wǎng)格可以較好地描述浮子的形狀。

表2 浮子周圍的網(wǎng)格數(shù)量和尺寸

圖3 比較了在不同網(wǎng)格尺寸下采用上述波浪條件捕捉浮子垂蕩運動的能力。結果表明，0.01 cm和0.005 cm 的網(wǎng)格尺寸下浮子位移相差不大，說明以0.01 cm 劃分的網(wǎng)格具有足夠高的模擬分辨率。因此，浮子周圍網(wǎng)格尺寸建議設置為Δx

圖3 不同網(wǎng)格尺寸的浮子垂蕩位移歷時曲線

2.3 三維數(shù)值水槽驗證

2.3.1 造波能力驗證

圖4 為距波浪入口2L附近（放置浮子處）規(guī)則波的數(shù)值計算結果與理論值的時程曲線（波高H=0.2 m，周期T=2 s，水深h=1 m）。在入射波傳遞一段時間到達穩(wěn)定狀態(tài)后，很明顯波峰處略尖陡，波谷處略平坦，且相對誤差低于4%；同時可以發(fā)現(xiàn)波浪的反射效應并不明顯，這表明水槽的造波能力和消波效果都比較可靠。

圖4 自由波面理論值與計算值的時程曲線

2.3.2 波浪對浮子作用的驗證

為驗證數(shù)值模型對浮子運動捕捉的準確性，選用Yu等[12]的試驗數(shù)據(jù)與數(shù)值計算結果進行比較。表3 給出了連續(xù)5 個振動周期中浮子運動最高點與最低點的數(shù)值結果與試驗數(shù)據(jù)的對比（波高H=0.2 m，周期T=2 s，水深h=1.0 m）?？梢钥闯?，數(shù)值結果與試驗數(shù)據(jù)的相對誤差均在10%以內(nèi)，所建立數(shù)值模型能夠較好地模擬波浪對浮子的作用。

表3 垂蕩位移數(shù)值結果與試驗數(shù)據(jù)

2.3.3 波浪沖擊壓力觀測點設置和定義

根據(jù)對稱性，在浮子迎浪方向、側面以及后方（分別為Row A、Row B 和Row C）設置3 列壓力傳感器，定義s為浪高儀至浮子底端的垂直距離，d為浮子高度，浮子柱面觀測點和壓力傳感器的布置如圖5（a）和表4 所示。浮子底部設置了11 個壓力傳感器，θ為所在點與x軸負方向的夾角，相鄰2 組測點之間的夾角為45°，浮子底面觀測點和壓力傳感器的布置如圖5（b）和表5 所示。模型中所有虛擬觀測點的布置與實際位置誤差小于3 mm。

表5 浮子底面壓力傳感器設置

圖5 波浪沖擊壓力虛擬觀測點設置

表4 浮子柱面壓力傳感器設置

波浪在沖擊的過程中涉及到諸多因素，比如水氣摻混、瞬時效應等，這些因素導致波浪沖擊垂蕩浮子的過程具有很強的隨機性，因此沖擊過程中的壓力峰值具有很強的離散性。本文根據(jù)文獻[18]對每個虛擬的壓力傳感器采集的壓力數(shù)據(jù)進行了處理，即浮子上的壓力測點取樣時，將每個周期內(nèi)的壓力峰值由大到小排列，并且取前1/3 個峰值的平均值作為各測點的沖擊壓力峰值特征值，即統(tǒng)計分析的特征值。每個測點的沖擊壓力峰值的實測數(shù)據(jù)記為Pt，浮子靜止時測點的初始壓力記為P0，將實測數(shù)據(jù)與初始壓力無量綱化后記為Pt/P0。浮子所受的水平和豎直方向波浪荷載峰值數(shù)據(jù)的提取則是將每個周期內(nèi)的壓力峰值提取出來，并且對峰值相差較小的幾組數(shù)據(jù)取平均值為波浪荷載的特征值。水平和垂直波浪荷載分別用Fx和Fz來表示，對波浪荷載無量綱化分析，選取浮子的重力ρgSd為無量綱化參數(shù)，其中ρ為浮子的密度，g為重力加速度，S為浮子底面積。

2.3.4 計算工況

數(shù)值計算中選用圓柱形浮子，實際浮子高度為1.6 m、質(zhì)量為3 860 kg、吃水深度為0.48 m。計算工況如表6 所示。

表6 計算工況設置

3 波浪對垂蕩浮子裝置沖擊作用

3.1 波浪沖擊壓力特征

3.1.1 沖擊壓力峰值在浮子的立側處分布規(guī)律

圖6 為工況1 中A3 傳感器記錄的浮子波浪壓力歷時曲線，由圖可以發(fā)現(xiàn)，浮子受到的沖擊壓力具有明顯的周期性，周期與入射波周期保持一致，沖擊壓力近似為正弦脈沖形式，沖擊階段壓力逐漸增大，之后隨著波面下降開始逐漸回落，當測點脫離水體，壓力維持在零值幾乎不變。

圖6 浮子立側不同壓力測點沖擊壓力歷時曲線

圖7 為工況1 條件下浮子表面不同測點上的壓力分布，1#和2#測點沖擊壓力在不同位置處差距較小，3#測點的沖擊壓力差別較大。因為波面高度達不到，故4#和5#測點的沖擊壓力為0。距離浮子底部越大，無量綱沖擊壓力值越大，即相較于初始壓力提升越大，3#測點的無量綱沖擊壓力達到最大值。距離浮子底部相同距離的測點中，Row C 側的沖擊壓力最大，Row A 側的沖擊壓力最小，表明浮子后側所受到的沖擊壓力最大。

圖7 浮子立側處不同測點壓力分布

3.1.2 沖擊壓力峰值在浮子底面分布規(guī)律

圖8 為工況1 條件下M1#、O#和S1#測點記錄的浮子的底部波浪壓力歷時曲線，由圖可以發(fā)現(xiàn)，浮子底端受到的沖擊壓力具有明顯的周期性，周期與入射波周期保持一致；浮子底部中心所受的波浪壓力小于波浪傳播方向上中心兩側的波浪壓力。

圖8 浮子底面不同壓力測點沖擊壓力歷時曲線

3.2 波浪參數(shù)和阻尼參數(shù)對浮子上波浪沖擊壓力垂直分布的影響

3.2.1 波陡參數(shù)對浮子周圍波浪沖擊效應影響

圖9 給出了不同波陡參數(shù)下，最大波浪沖擊壓力沿的垂直分布。

從圖9 中可以發(fā)現(xiàn)，最大波浪沖擊壓力隨著波陡參數(shù)的增加而增大，距離浮子底部距離越近的位置上波浪沖擊壓力越大。波陡較大時，浮子靠上位置處的波浪沖擊壓力最小值出現(xiàn)在Row B 側，這是由于浮子側面的波浪爬升高度處于最低導致，其他位置的最大沖擊壓力出現(xiàn)在浮子中間（Row C），最小沖擊壓力依然發(fā)生在迎浪側（Row A）。

圖9 不同波陡參數(shù)下浮子上最大壓力垂直分布

3.2.2 散射參數(shù)對浮子周圍波浪沖擊效應影響

圖10 給出了不同散射參數(shù)下，最大波浪沖擊壓力沿Row A、Row B 和Row C 的垂直分布。從圖10 可以看出，改變浮子半徑對于Row B 和Row C側相同壓力測點上波浪沖擊壓力峰值的改變較小。散射參數(shù)較小時，Row A 側的相同測點壓力峰值隨散射參數(shù)變化較??；散射參數(shù)較大時，Row A側的相同測點壓力峰值有明顯提升?？偟膩碚f，增加浮子半徑對于浮子截面上相同位置的波浪沖擊壓力影響較小。

圖10 不同散射參數(shù)下浮子上最大壓力垂直分布

3.2.3 線性PTO 阻尼參數(shù)對浮子周圍波浪沖擊效應影響

圖11 給出了不同彈簧剛度系數(shù)下，最大波浪沖擊壓力沿Row A、Row B 和Row C 的垂直分布。

從圖11 中可以看出，彈簧剛度系數(shù)的增加對于浮子上的波浪沖擊壓力影響較小，浮子上相同壓力測點上的壓力峰值隨彈簧系數(shù)的增加幾乎保持不變。通過對比發(fā)現(xiàn)，施加PTO 后對于浮子上的波浪沖擊壓力改變很小，且浮子上的波浪沖擊壓力垂直分布規(guī)律與浮子自由垂蕩運動時保持一致。

圖11 不同彈簧剛度系數(shù)下浮子上最大壓力垂直分布

圖12 給出了不同阻尼系數(shù)下，最大波浪沖擊壓力沿Row A、Row B 和Row C 的垂直分布。從圖12 中可以看出，增加阻尼系數(shù)對浮子截面上波浪沖擊壓力垂直分布的影響較??；增加PTO 阻尼對浮子垂蕩運動有影響，但浮子上的波浪沖擊由于波浪參數(shù)保持一致，導致浮子上的波浪沖擊壓力受其他因素的影響很小。

圖12 不同阻尼系數(shù)下浮子上最大壓力垂直分布

浮子側面的波浪沖擊壓力受波陡參數(shù)的影響最大，沖擊壓力隨波陡增大而增大；Row A 側的波浪沖擊壓力隨散射參數(shù)增大而增加，其他位置處的沖擊壓力受散射參數(shù)的影響不明顯；浮子上的最大波浪沖擊壓力受PTO 系統(tǒng)阻尼的影響很??；浮子側面的壓力分布在s>0.24 m 位置上的壓力處于最小值，這主要是由于波浪在該位置的爬升較小，沖擊強度較低所導致。

3.3 最大波浪荷載與波浪沖擊壓力關系研究

圖13 和圖14 給出了工況3 水平和垂直波浪荷載與Row A 壓力測點時程對照。A1 和A2 處的波浪壓力與水平波浪荷載的歷時曲線存在較小的相位差，水平波浪荷載的最大值與最小值出現(xiàn)的時間較早于A1 處的波浪壓力最值；浮子迎浪側距離底部較高位置處（A3、A4 和A5）的波浪壓力與水平波浪荷載的歷時曲線幾乎不存在相位差，即水平波浪荷載的最大值與最小值出現(xiàn)的時間和該位置處的波浪沖擊壓力最值同步出現(xiàn)。

圖13 水平波浪荷載與Row A 壓力測點時程曲線對照

圖14 垂直波浪荷載與Row A 壓力測點時程曲線對照

由于水平波浪荷載受力位置與Row A 側相似，都為迎浪面，因此歷時曲線存在的相位差較小，在浮子較高位置處的相位差幾乎一致。但浮子所受垂直波浪荷載的受力面與浮子Row A 側波浪壓力受力面存在較大的差異，從圖14 中可以發(fā)現(xiàn)，垂直波浪荷載的最值出現(xiàn)時間和該位置處的波浪沖擊壓力最值出現(xiàn)時間存在較大差異。A1 和A2 處的波浪壓力與垂直波浪荷載最值出現(xiàn)的時間相反，即垂直波浪荷載出現(xiàn)最大值時，該位置的波浪壓力處于最小位置；浮子迎浪側距離底部較高位置處（A3、A4 和A5）的波浪壓力與垂直波浪荷載的最值同步出現(xiàn)，與水平荷載相似。

4 結論

1）在規(guī)則波作用下，浮子所受沖擊壓力具有明顯的周期性，浮子底部波浪沖擊壓力近似為正弦脈沖形式。

2）在入射波前進方向上，隨著浮子底部傳感器與浮子前部距離的增大，傳感器所受的沖擊壓力峰值也增大。

3）波浪沖擊壓力受波陡參數(shù)的影響最大，沖擊壓力隨波陡增大而增大。浮子迎浪側的波浪沖擊壓力隨散射參數(shù)的增大有增加的趨勢，其他位置處的沖擊壓力受散射參數(shù)的影響不明顯。浮子上的最大波浪沖擊壓力受PTO 阻尼系統(tǒng)的影響很小。

4）浮子迎浪側距離浮子底部較高位置處的波浪沖擊壓力與波浪荷載最值出現(xiàn)時間幾乎保持一致。

5）本研究只針對線性波浪作用下的浮子沖擊壓力進行了研究，隨機波浪作用下浮子荷載特征的研究需進一步開展。