張 亮，國 威，王樹齊

(哈爾濱工程大學船舶工程學院，150001哈爾濱)

波浪能作為一種清潔的可再生能源，蘊藏量巨大，為了緩解能源危機，波浪能裝置的研發(fā)引發(fā)人們的日益關(guān)注.波能裝置按結(jié)構(gòu)形式分為點頭鴨式［1］、振蕩水柱式［2-3］、推擺式［4-5］、聚波蓄能式［6］、點吸式［7-9］.點吸式波能裝置有建造難度低、成本低、效率高等優(yōu)點，已成為國內(nèi)外學者研究的熱點.平麗［10］基于線性勢流理論，用邊界元法研究裝置的水動力性能，得出裝置的俘獲寬度比和最優(yōu)阻尼的表達式，該方法忽略了流體的黏性效應.歷福偉［11］利用ADINA軟件建立了二維數(shù)值波浪水池，研究了浮子形狀對振蕩浮子式發(fā)電裝置效率的影響.該方法的缺點是忽略了流體黏性，二維模擬僅限于形狀簡單的浮子，適用范圍較窄.Guanche等［12］利用Simulink對一種波能裝置進行了部分非線性時域數(shù)值模擬，通過改變電機的阻尼系數(shù)、裝置直徑、水深、軸承摩擦系數(shù)、垂向拖曳力系數(shù)、波高、周期等研究裝置的水動力性能，此模型形狀較為簡單，沒有對形狀略復雜的浮子進行分析，不利于不同形狀浮子水動力性能的對比.Vicente［13］對一種帶有張緊錨鏈的點吸式波能裝置的非線性水動力進行了時域和頻域分析.Bhatta［14-15］等運用勢流理論對處于有限水深的圓柱體的繞射和輻射問題進行了數(shù)學推導，并得出數(shù)值解，得到了圓柱體的在規(guī)則波中的運動規(guī)律.目前對點吸式波浪能裝置的理論研究多以勢流理論為主，裝置中浮子形狀過于簡單，而且忽略了流體的黏性效應，與實際情況存在差異.

本文運用CFD方法，針對一種特定的點吸式波能裝置進行三維數(shù)值模擬，對浮子形狀及外界工作參數(shù)進行了優(yōu)化分析.本方法的優(yōu)勢在于考慮了流體的黏性，可以對形狀較為復雜的浮子進行計算，比勢流理論更接近真實情況.

1 理論分析

1.1 浮子運動的數(shù)學模型

點吸式波能裝置可簡化為質(zhì)量-彈簧-阻尼系統(tǒng)，浮子在波浪中垂蕩運動時，受到慣性力、阻尼力、回復力和波浪力的作用，其中阻尼力與浮子垂蕩速度成正比，回復力與浮子位移成正比.建立浮子的垂向受力平衡關(guān)系式:

式中:2ν為垂蕩運動阻尼系數(shù)，2ν=(C+C0)/(m+λ);n為浮子垂蕩固有頻率，n2=(ρgs+K)/(m+λ);A為波幅;X0為波幅修正系數(shù);ω為波浪頻率;C為機械阻尼系數(shù);C0為水的垂蕩阻尼系數(shù);m為浮子質(zhì)量;K為彈性系數(shù).

浮子垂蕩的固有周期

由式(1)得出浮子運動幅值z與波幅A的關(guān)系如下:

式中:Λ為波浪頻率與浮子固有頻率之比，Λ=ω/n;μ為無因次衰減系數(shù).

入射波的功率Pi為浮子直徑寬度內(nèi)的波浪功率，且

式中:H為波高，T為波浪周期，D為浮子直徑.浮子的平均吸收功率為

式中:FG為電機作用在浮子上的阻尼力，v為浮子垂蕩速度，t2-t1=nT，n為整數(shù).

浮子的俘獲寬度比η定義為浮子平均輸出功率與浮子寬度內(nèi)波浪輸入功率之比，即

浮子在做垂蕩運動時對電機做功，通過電機的阻尼力轉(zhuǎn)化為電能.

1.2 CFD數(shù)學模型

CFD方法基于黏性流體力學理論，采用基于有限元的有限體積法，在保證了有限體積法的守恒特性的基礎(chǔ)上，吸收了有限元法的數(shù)值精確性，如式(6)和式(7)所示.

連續(xù)性方程:

RANS方程組:

式中:ui和uj代表速度分量時均值，i，j=1，2，3;p為壓力時均值;ρ為流體密度;ν為流體的運動黏性系數(shù);gi為重力加速度分量為雷諾應力項.

針對本文的裝置模型，選取合適的邊界條件，采用動網(wǎng)格技術(shù)模擬浮體運動，采用體積分數(shù)法標記自由面，選取k-ω湍流模型，該湍流模型具有計算準確，收斂性好的特點，如式(8)和式(9):

k方程為

ω方程為

式中:k為湍動能，ω為湍流頻率，Pk為湍流生成速率，β'=0.09，α=5/9，β=0.075.

2 數(shù)值計算

2.1 計算模型

運用ANSYS13.0中的CFX軟件對浮子運動進行模擬，計算模型采用實體尺寸，如圖1(a)所示，浮子上部為圓柱，下部為圓臺，中間開孔.計算域的大小為50 m×20 m×50 m，采用推板造波，通過動網(wǎng)格模擬造波板運動，生成規(guī)則波［16］，造波板、水池池壁、底部設置為No Slip Wall，水池頂部設置為Opening，水池尾部為消波區(qū)，可以大幅減少波浪反射.

計算域網(wǎng)格的劃分如圖1(b)所示，在波長方向?qū)W(wǎng)格數(shù)量的控制取決于單位波長內(nèi)的網(wǎng)格的數(shù)量，實踐證明一個波長內(nèi)應有120～150個網(wǎng)格，波高方向上30～50個網(wǎng)格，既保證了數(shù)值計算的精準度，又節(jié)省了計算時間.

圖1 浮子尺寸及計算域網(wǎng)格示意

浮子位置距造波板15 m，浮子的運動求解采用剛體求解，約束其他自由度，使浮子僅做單自由度的垂蕩運動，浮子力和運動的收斂標準取為10-4，計算時每個時間步的迭代次數(shù)在10～15次.時間步長為0.05 s，模擬10～15個波浪周期下浮子的運動.

2.2 有效性驗證

為了驗證CFD方法的有效性與準確性，分別采用商業(yè)CFD軟件CFX和基于勢流理論的計算軟件AQWA對實驗模型進行數(shù)值模擬，并與實驗值對比，如圖2所示.

實驗模型的縮尺比為10，質(zhì)量為14.25 kg，該實驗在哈爾濱工程大學拖曳水池完成，實驗波高為 0.08 m，周期為 2.2 s.

從圖2中可以看出，運用CFD方法對浮子垂蕩位移和速度的模擬結(jié)果與實驗結(jié)果吻合良好，驗證了CFD方法的有效性;將CFD值與勢流理論值進行對比表明，CFD值比勢流理論值更接近實驗值，驗證了該方法的準確性.綜上可知，CFD方法對浮子運動的模擬有效且準確.

3 結(jié)果分析

波浪能裝置的性能與很多因素相關(guān)，模擬浮子在不同形狀、不同電機阻尼、不同質(zhì)量下的運動，并分析得出了浮子在不同外界參數(shù)下的性能.

3.1 圓柱和圓臺浮子的性能對比

計算分析兩種形狀的浮子在不用波況下的水動力性能，如圖3所示，分為浮子在無負載和有負載下兩種工況.

圖3 圓臺和圓柱形浮子形狀示意

3.1.1 自然振蕩

無負載時浮子在在波浪中做垂蕩運動，此時浮子所受阻尼力來自于流體，裝置阻尼系數(shù)僅為C0.式(2)說明了浮子的運動幅值與波幅的關(guān)系，將式(2)變形得

式中:K為放大系數(shù)，修正系數(shù)X0為1.

選取質(zhì)量為17.528 t，直徑和高度相同的圓柱和圓臺進行對比，波高為0.4 m，波浪周期在2.40～5.06 s，如圖 4 所示.

圖4 圓臺和圓柱在不同波浪周期下無機械阻尼時的幅值響應

從圖4可以看出，隨著波浪周期的增大，兩者的運動幅值響應先增大后減小，原因如下:

根據(jù)式(10)可知，當浮子的固有周期大于波浪周期，即Λ＞1時，浮子的運動受波浪影響很小，當波浪周期逐漸增大到Λ=1時，K出現(xiàn)最大值，此時浮子固有周期與波浪周期相同，發(fā)生共振，當波浪周期大于浮子固有周期即Λ＜1時，浮子隨波浪漂浮，波浪周期足夠大時，浮子的運動幅值與波浪幅值相同，K趨近于1.

由于兩種浮子的形狀不同，其共振周期和K值都略有不同，圓臺的共振周期為3.7 s，圓柱的共振周期為3.3 s，除圓柱的共振點附近的區(qū)域外，在波浪周期2.40～5.06 s范圍內(nèi)，圓臺的運動幅值都高于圓柱，這說明同尺度、同質(zhì)量的圓柱和圓臺，圓臺對波浪的響應要優(yōu)于圓柱.

3.1.2 阻尼振蕩

上述分析是圓柱和圓臺浮子在無負載下做垂蕩運動，裝置在有外界機械阻尼時，阻尼系數(shù)增大，使無因次衰減系數(shù)μ增大，根據(jù)式(10)可知，浮子運動幅值將減少.圖5為圓柱和圓臺在機械阻尼為10 000 Ns/m時的幅值響應曲線.

圖5 圓臺和圓柱在不同波浪周期下有機械阻尼時的幅值響應

從圖5中可以看出，在相同波況下圓柱的運動幅值小于圓臺，這是由于在波況相同時，圓臺垂蕩時作用于周圍流體而產(chǎn)生的能量耗散小于圓柱，圓臺自身剩余能量較大，其運動幅值也較大.從圖中還可以得出圓柱的共振周期為3.5 s，而無阻尼時其共振周期為3.3 s，兩者略有差異，這是由于浮子的共振與固有周期有關(guān)，固有周期的大小與浮子的附加質(zhì)量有關(guān)，由于浮子在有阻尼和無阻尼時運動幅值不同，其周圍流場存在差異，導致附加質(zhì)量會有變化，從而影響其固有周期.

根據(jù)上述分析可知，在自然振蕩和阻尼振蕩時，圓臺形浮子的垂蕩運動性能要優(yōu)于圓柱形浮子.

3.2 機械阻尼對浮子俘獲寬度比的影響

圓臺浮子運動性能較圓柱好，對圓臺浮子進行模擬，選取波高為0.4 m，圓臺浮子的共振周期為3.7 s的波況，浮子的質(zhì)量和尺寸不變，對浮子在機械阻尼為10 000～200 000 Ns/m進行計算，圖6(a)給出了浮子在機械阻尼系數(shù)為10 000、30 000和50 000 Ns/m下的運動曲線.

圖6 不同阻尼系數(shù)下的浮子垂蕩速度和俘獲寬度比曲線

從圖6(a)可以看出，在給定波況下，浮子的垂蕩速度曲線趨于穩(wěn)定后，其運動幅值基本保持不變，機械阻尼系數(shù)越大，浮子的速度幅值越小.通過分析可知，電機的阻尼系數(shù)為零時，浮子處于自由振蕩狀態(tài)，所受電機阻尼力為零，理論上浮子的輸出功率為零;當電機阻尼系數(shù)無限大時，相當于浮子被剛性固定，速度為零，其輸出功率也為零，因此找到最佳電機阻尼是必要的.圖6(b)為不同阻尼系數(shù)下浮子俘獲寬度比曲線，從圖中可以看出，浮子的俘獲寬度比隨著機械阻尼的增大先增大后減小，在阻尼系數(shù)為70 000 Ns/m時，俘獲寬度比達到最大值0.8.曲線先上升后下降的原因在于:根據(jù)式(3)～(5)可知，浮子寬度內(nèi)的波浪功率與波高、周期和浮子直徑有關(guān)，當浮子處于給定的波況時，波浪功率為定值，浮子俘獲寬度比只與電機阻尼力和速度大小有關(guān)，浮子瞬時吸收功率為

由式(11)可知，浮子的瞬時吸收功率大小與阻尼系數(shù)和速度的平方成線性關(guān)系，而浮子的速度大小受電機阻尼系數(shù)的影響，浮子的運動速度隨阻尼的增加而減少，阻尼越大，其克服外界阻尼力做功越多，導致其自身動能減少，速度就越小.總體表現(xiàn)為輸出功率在某個速度上出現(xiàn)最大值.

3.3 質(zhì)量對浮子俘獲寬度比的影響

在對浮子形狀和機械阻尼進行優(yōu)化后，進一步計算圓臺浮子質(zhì)量對俘獲寬度比的影響.選取波高為 0.4 m，周期為 4.4 s，機械阻尼系數(shù)為30 000 Ns/m，對不同質(zhì)量的浮子進行計算.選取浮子質(zhì)量范圍為11～18 t，計算得出浮子質(zhì)量-俘獲寬度比曲線，如圖7所示.

圖7 浮子的質(zhì)量-俘獲寬度比曲線

圖7中不同浮子質(zhì)量對應不同的吃水，隨著浮子質(zhì)量的增加，其俘獲寬度比先增大后緩慢減小，這是由于浮子質(zhì)量增加后其吃水增加，使浮子所受波浪力對浮子做功增加，但隨著吃水的繼續(xù)增大會使浮子的慣性也相應增大，浮子對波浪的響應度降低，導致俘獲寬度比略有下降.但在不同浮子質(zhì)量的情況下，該波能裝置的俘獲寬度比浮動不超過2%.

由此可見，在遠離浮子共振周期時，浮子質(zhì)量的改變對裝置的性能影響很小.在滿足性能要求的基礎(chǔ)上可適當減少浮子質(zhì)量.

4 結(jié)論

1)相同尺度和質(zhì)量的圓柱和圓臺浮子的共振周期略有不同，自然振蕩下的浮子的運動幅值要高于阻尼振蕩時的運動;在兩種振蕩條件下，圓臺浮子的運動性能優(yōu)于圓柱浮子.

2)隨著機械阻尼的增大，裝置的俘獲寬度比先增大后減小，表明裝置存在最佳的機械阻尼，機械阻尼過大或過小都不利于浮子對能量的吸收.

3)質(zhì)量的變化對俘獲寬度比的影響不明顯，考慮制造成本，可適當調(diào)節(jié)浮子質(zhì)量.

4)基于CFD方法對浮子運動進行三維數(shù)值模擬，并對浮子在不同參數(shù)下的性能進行分析，得出最佳工況，實現(xiàn)了對該特定點吸式波能裝置性能的優(yōu)化.

［1］余志.海洋波浪能發(fā)電技術(shù)進展［J］.海洋工程，1993，11(1):86-93.

［2］韓冰峰，褚金奎，熊葉勝，等.海洋波浪能研究進展［J］.電網(wǎng)與清潔能源，2012，28(2):61-66.

［3］CLéMENT A，MCCULLEN P，F(xiàn)ALCO A.Wave energy in Europe:currentstatus and perpectives ［J］.Renewable＆ Sustainable Energy Reviews，2002(6):405-431.

［4］李允武.海洋能源開發(fā)［M］.北京:海洋出版社，2008.

［5］CAMERON L，DOHERTY R，HENRY A，et al.Design of the next generation of the Oysterwave energy converter［C］//3rd International Conference on Ocean Energy.Bilbao:Amanda W，2010:22-39.

［6］任建莉，鐘英杰，張雪梅，等.海洋波能發(fā)電的現(xiàn)狀與前景［J］.浙江工業(yè)大學學報，2006，34(1):69-73.

［7］ZURKINDEN A S，F(xiàn)ERRI F，BEATTY S，et al.Nonlinear numerical modeling and experimental testing of a pointabsorberwave energy converter［J］.Ocean Engineering，2014，78(1):11-21.

［8］SOREN R K，BISWAJIT B，MAHDI T，et al.Optimal control of an array of non-linear wave energy point converters［J］.Ocean Engineering，2014，88(15):242-254.

［9］LOPES M F P，HALS J，MOAN T，et al.Experimental and numerical investigation of non-predictive phase-control strategies for a point-absorbing wave energy converter［J］.Ocean Engineering，2009，36(5):386-402.

［10］平麗.振蕩浮子式波能轉(zhuǎn)換裝置性能的研究［D］.大連:大連理工大學，2005.

［11］歷福偉.浮子形狀對振蕩浮子式波浪發(fā)電裝置效率的影響［D］.哈爾濱:哈爾濱工業(yè)大學，2011.

［12］GUANCHE R，GOMEZ V，VIDAL C，et al.Numerical analysis and performance optimization of a submerged wave energy point absorber［J］.Ocean Engineering，2013，59(1):214-230.

［13］VICENTE P C.Nonlinear dynamics of a tightly moored point-absorberwave energy converter［J］.Ocean Engineering，2013，59(1):20-36.

［14］BHATTA D D，RAHMAN M.On scattering and radiation problem for a cylinder in water of finite depth［J］.International Journal of Engineering Science，2003，41(9):931-967.

［15］ BHATTA D D. Computation ofhydrodynamic coefficients，displacement-amplitude ratios and forces for a floating cylinder due to wave diffraction and radiation［J］.International Journal of Non-Linear Mechanics，2011，46(1):1027-1041.

［16］王永學.無反射造波數(shù)值波浪水槽［J］.水動力學研究與進展，1994，9(2):205-214.