孫科， 解光慈， 周斌珍

(1.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001; 2.中山大學(xué) 海洋工程與技術(shù)學(xué)院，廣東 珠海 519082； 3.華南理工大學(xué) 土木與交通學(xué)院,廣東 廣州 510000)

波浪能作為一種新型的海洋可再生能源，具有儲量大、分布廣的特點[1]。而波浪能裝置的形式也是多種多樣的，有點吸收式、振蕩水柱式以及越浪式等[2]。如今，越來越多的國家研究波浪能利用技術(shù)，主要研究方面包括波能轉(zhuǎn)換水動力研究、PTO系統(tǒng)研究以及波能浮子與平臺的耦合運動等[3]。其中浮子的形狀很大程度上決定了其水動力性能。

對于波浪能浮子的研究方法有很多種，最常用的方法有基于勢流理論的邊界元法和基于N-S方程的計算流體動力學(xué)(CFD)方法。

在波浪的研究中，勢流理論一直占據(jù)主導(dǎo)地位。2013年，周斌珍等[4]利用高階邊界元研究波浪與結(jié)構(gòu)物的非線性相互作用；張萬超等[5]在2016年基于勢流理論對軸對稱浮子進行水動力特性研究。Zhang等[6]基于勢流理論和高階邊界元法，建立了只考慮雙浮體垂蕩運動的數(shù)值模型，研究其水動力性能和時域響應(yīng)。常見的基于勢流理論的軟件有WAMIT、AQWA以及SEASUM。然而由于其無旋無粘的前提假設(shè)，使得這一理論只適用于處理一些大結(jié)構(gòu)物或者粘性影響并不顯著的流動問題。

而隨著計算機性能的不斷發(fā)展，計算流體動力學(xué)隨之誕生，通過對整個計算域離散化，利用有限差分法、有限體積法等方法計算。CFD方法能夠考慮物體的粘性，并且可對復(fù)雜幾何形狀的物體計算，因此成為許多研究者們的選擇。國威[7]在2016年基于CFD方法模擬浮子的線性波上浪、非線性波抨擊及自由面形變的時域非線性問題；Bhinde等[8]在2009年基于CFD方法研究一種點吸收式的波能浮子的水動力特性以及幾何形狀的優(yōu)化并與實驗作對比；Ransley等[9]在2017年基于CFD方法研究波能浮子在極端海況的生存能力；Wang[10]通過綜述近年來用CFD數(shù)值方法解決的船舶與海洋工程中的復(fù)雜問題，討論了立管渦激振動、浮式風(fēng)力機以及高速船在波浪中的運動等方面的引用，來描述計算流體力學(xué)技術(shù)的發(fā)展。常見的CFD商用軟件有CFX、STARCCM+以及FLUENT等。然而CFD方法對于網(wǎng)格的數(shù)量及質(zhì)量要求很高，計算時間長，對計算機性能要求高，不適合大范圍的算例計算。

對于波浪能裝換裝置來說，浮子選型對能量轉(zhuǎn)換效率有很大影響，國內(nèi)外學(xué)者針對這一方面做出了大量研究。張亮等[11]基于CFD方法研究圓柱和圓臺浮子的阻尼系數(shù)、浮子質(zhì)量對裝置性能的影響規(guī)律。Zhang等[12]通過邊界離散化(BDM)的半解析方法分析復(fù)雜軸對稱浮體，分析對比了圓柱、錐形以及半球形浮子的水動力性能。然而以上研究對于浮子形狀分類較少以及沒有結(jié)合分析流場機理。而關(guān)于粘性修正方面，同樣也有許多學(xué)者在研究，Liu等[13]在對浮體進行分析時，通過在自由表面邊界中引入耗散項來考慮粘性效應(yīng)。田博等[14]基于RANS所求得的水動力系數(shù)與勢流理論結(jié)合進行粘性修正，對單體復(fù)合船型運動預(yù)報方法進行研究。

本文利用CFD與勢流理論結(jié)合的方法，能夠減少大量計算時間。通過CFD模擬自由衰減估算粘性力，再以勢流理論的方法加上粘性力對波能浮子的幾何形狀進行研究，分析平底圓柱浮子的直徑吃水比對粘性力、運動響應(yīng)及俘獲寬度比的影響；對比不同底部角度錐形浮子運動響應(yīng)及俘獲寬度比，選擇最優(yōu)寬角度。

1 波能浮子數(shù)學(xué)模型的建立

當(dāng)波浪傳遞時，波能浮子受到波浪作用產(chǎn)生垂蕩作用，使得浮子與PTO產(chǎn)生相對的位移，進而產(chǎn)生能量。因此在數(shù)學(xué)模型建立時只考慮浮子垂蕩方向的運動。對于圓柱形波能浮子，質(zhì)量為m，半徑為a,如圖1。

圖1 波能浮子運動模型

對于單自由度波能浮子，可用水面浮體的運動表示，建立運動方程為：

(1)

(2)

Fs=Cx=ρgπa2x

(3)

(4)

Fk=Kx3

(5)

式中：Fe波浪激勵力；Fr輻射力；Fs靜水回復(fù)力；Fc線性PTO阻尼力；Fk彈簧力；μ33浮子附加質(zhì)量；λ33浮子輻射阻尼；ρ為水的密度；g為重力加速度；CpPTO阻尼力系數(shù)；K為彈簧系數(shù)。

將各項代入得：

(6)

在線性勢流理論中，浮子在波浪中運動時所受的力可分為入射力、繞射力以及輻射力、入射力和繞射力統(tǒng)稱為波浪激勵力。在實際水域中，浮子由于受到粘性影響，其運動響應(yīng)尤其在共振頻率下遠遠小于勢流理論所得結(jié)果，流體粘性力是導(dǎo)致這一結(jié)果的主要原因。Son等[15]對圓柱浮子進行波浪能裝置實驗時發(fā)現(xiàn)，勢流理論所得波浪激勵力結(jié)果與實驗所測基本相同，因此，波浪激勵力項受粘性影響不大。但對于輻射力項，尤其是阻尼項，試驗和勢流的結(jié)果相差很大，粘性效應(yīng)不可忽略。為了能夠得到浮子的粘性力，Journée等[16]曾通過自由衰減曲線推導(dǎo)浮子的粘性力公式，當(dāng)浮子僅考慮垂蕩運動時，其垂蕩自由衰減線性運動方程為：

(7)

化簡后得：

(8)

設(shè)無量綱系數(shù)k：

(9)

(10)

(11)

x3(t)=x3,ae-vt[cos(ωvt)+sin(ωvt)]

(12)

設(shè)x3(t)周期為Tv，求解：

x3(t)/x3(t+Tv)=evTv

(13)

化簡得：

vTv=kω0Tv=ln(x3(t)/x3(t+Tv))

(14)

ω0Tv=ωvTv

(15)

繼而可以推導(dǎo)無量綱阻尼：

(16)

為了防止由于零點漂移原因使得所得數(shù)據(jù)不準，通過計算衰減曲線的雙振幅數(shù)值，來避免這種誤差。

(17)

式中x3,k為第k個峰值點的幅值。

式(17)能夠通過自由衰減曲線得到無量綱系數(shù)k，進而估算粘性阻尼系數(shù):

λvis=2kρgπa2/ω0

(18)

最終推導(dǎo)出浮子固有頻率時的粘性阻尼，Tom[17]在實驗中證明，波浪頻率變化對粘性效應(yīng)影響不大，因此可以通過以上方法算出固有頻率下的總阻尼系數(shù)。結(jié)合Tom的結(jié)論，可以大致計算出粘性阻尼，從而對整個頻域進行修正。

浮子粘性阻尼計算：

λv=λvis-λn

(19)

因此，通過上述論述，得到浮子在垂直方向的運動方程：

(20)

式中λn為共振頻率下浮子輻射阻尼。

每個波浪頻率下，浮子的輸出功率存在最優(yōu)阻尼系數(shù)，即cp為某一值時，輸出功率最大：

(21)

2 數(shù)值模擬方法

利用CFD軟件Star-CCM+來模擬浮子自由衰減，進而算出粘性系數(shù)。模型設(shè)置的是一個四周壁面的三維水池，浮子自由衰減所興的波會向四周傳遞，由于水池大小有限，需要考慮壁面效應(yīng)，因此設(shè)置阻尼消波。

模擬中涉及到空氣與水的多項流，應(yīng)用流體體積(VOF)法捕捉空氣相和水相之間的自由表面，因此需要對液面附近做網(wǎng)格加密，浮子運動區(qū)域也同樣需要設(shè)置網(wǎng)格加密，并且為了網(wǎng)格尺寸能夠平滑過渡，設(shè)置過渡層。網(wǎng)格設(shè)置如圖2所示：

圖2 CFD模型網(wǎng)格劃分

網(wǎng)格劃分及模型設(shè)置是CFD前處理必不可少的，理論上來講，網(wǎng)格密度越高，計算精度越高，但同樣計算時間會增加。時間步長對于模擬是否能夠收斂也有很大影響。合適的網(wǎng)格尺寸和時間步長對計算精度和計算效率改善有著很好的效果。本文設(shè)置了3個測試模型，選擇2a/d=1.6的圓柱浮子，d為浮子吃水，T為浮子共振周期。其時間步長和網(wǎng)格大小設(shè)置如表1所示，模擬結(jié)果如圖3所示，根據(jù)3個模型數(shù)據(jù)得出粘性阻尼，并與模型1對比得相對誤差，如表2所示。由表2可知當(dāng)網(wǎng)格及時間步長選擇d/30和T/400時，所得結(jié)果與另外2組相差不大，數(shù)據(jù)差異小于5%，因此d/30和T/400的網(wǎng)格及時間步長滿足收斂性要求。

表1 自由衰減模型時間步長與網(wǎng)格尺寸設(shè)置

自由衰減模擬除了對網(wǎng)格大小及時間步長收斂性驗證以外，還要對其進行準確性驗證。Tom[17]研究了圓柱浮子的水動力性能，并做了相關(guān)浮子自由衰減實驗，將CFD數(shù)值模擬所得數(shù)據(jù)與其試驗數(shù)據(jù)作對比。

圖3 網(wǎng)格尺寸和時間步長收斂性驗證

表2 粘性阻尼相對誤差

觀察圖4，模擬數(shù)據(jù)與實驗數(shù)據(jù)，在周期上幾乎吻合，而在衰減曲線上，實驗數(shù)據(jù)后幾個周期略大于模擬數(shù)據(jù)，數(shù)據(jù)差別小于5%，驗證了數(shù)值模型的準確性。

圖4 自由衰減運動驗證

根據(jù)得到的粘性系數(shù)以及推導(dǎo)的運動方程，求解浮子的頻域運動響應(yīng)。為了驗證其修正數(shù)值的精度，將其與CFD模擬的波浪下浮子運動作對比，模型選擇2a/d=1.6的圓柱浮子，周期T=2.17 s，波高H=0.2 m，所設(shè)置阻尼均為最優(yōu)阻尼，得到圖5。

圖5 CFD模擬與粘性修正方法對比

綜合以上數(shù)據(jù)對比，基于粘性修正的勢流理論方法與CFD模擬相比，在共振頻率附近有輕微的差別，垂蕩運動響應(yīng)結(jié)果基本一致，說明可以通過粘性修正方法來估算浮子運動，減少計算量，增加計算效率。

3 平底圓柱浮子的水動力性能分析

由于平底圓柱形浮子幾何特征只有直徑2a和吃水d，因此幾何特征形狀用無量綱參數(shù)2a/d來表示。通過改變直徑吃水比，來改變浮子幾何尺度，更加全面地分析不同尺度浮子的水動力性能。為了更好地觀察不同直徑吃水比下浮子的粘性影響，定義阻尼粘性效應(yīng)的無量綱表示為：

(22)

通過STARCCM++模擬浮子的自由衰減，代入粘性修正公式，并將其無量綱化，得到不同2a/d下浮子的粘性系數(shù)變化。

圖6 平底圓柱浮子阻尼修正系數(shù)隨2 a/d變化關(guān)系

圖7中浮子的運動響應(yīng)峰值隨2a/d的增加而減小，當(dāng)頻率過大或者過小時，不同比例的浮子運動響應(yīng)差距不大。當(dāng)ω>2 rad/s左右時，運動響應(yīng)開始上升，且2a/d越大上升得越早，同時下降得也早；當(dāng)ω>5 rad/s左右時，所有比例基本重合。

圖8中浮子俘獲寬度比(CWR)隨著2a/d的增大其峰值逐漸增加，其所俘獲大功率的頻率范圍也更大，即有效頻域?qū)挾雀?。此外隨著2a/d的增加，浮子的峰值所處頻率也在不斷減小，當(dāng)2a/d=0.4時，峰值點處ω=2.9 rad/s，當(dāng)2a/d=3.6時，峰值點處ω=2 rad/s，可以通過調(diào)節(jié)尺寸選擇適合海域的浮子。大部分海域都有一段頻繁出現(xiàn)波浪周期。根據(jù)多發(fā)的波浪頻率設(shè)計相應(yīng)尺寸的浮子，對于工程應(yīng)用具有重大意義。

圖7 圓柱浮子頻域運動響應(yīng)

圖8 圓柱浮子頻域俘獲寬度比

4 錐底浮子的水動能性能分析

分析完浮子2a/d對平底浮子水動力性能的影響，接下來分析錐底浮子底部形狀對水動力性能的影響，選擇4種浮子進行分析，分別是α=90°,120°,150°,180°(即平底圓柱浮子),如圖9所示。

圖9 浮子形狀

選擇尺寸2a/d=1.6，4個浮子排水體積相同，錐底浮子的等效吃水與平底圓柱浮子的相同。

由圖10可知阻尼修正系數(shù)隨著角度的增大而增大，即底部形狀的角度越小所受粘性影響越小。計算浮子的運動響應(yīng)以及俘獲寬度比，由圖11可以觀察到，不同錐底角度的浮子，其共振頻率有細微的變化，當(dāng)角度變小時，因為小角度錐底浮子所受粘性影響小，其運動響應(yīng)就隨之增大。觀察俘獲寬度比變化圖可知隨著角度的減小，其俘獲寬度比的峰值增大，且有效頻域?qū)挾纫搽S之增大，更加有利于對波浪功率的俘獲,如圖12所示。然而在工程實際利用中，當(dāng)錐底角度越來越小，為保證其排水體積不變，長度就會隨之增加，不利于浮子的安裝以及整個裝置的運行，因此，對于適合角度的錐底浮子需要根據(jù)實際情況決定。

圖10 浮子阻尼修正系數(shù)隨底部角度變化關(guān)系

圖11 不同底部形狀浮子運動響應(yīng)

圖12 不同底部形狀浮子俘獲寬度比

在使用粘性修正的方法計算浮子水動力系數(shù)后，為保證數(shù)值模擬的準確性，需要利用CFD方法進一步分析計算，通過模擬浮子在波浪下的運動流場，分析流場結(jié)構(gòu)和變化機理，進一步檢驗線性計算結(jié)果的可靠性。本文利用Star-CCM+模擬4個浮子在共振周期附近的波浪下的運動，監(jiān)測其流場渦的變化，波浪為從右向左傳遞。

由圖13～16為4個浮子在在一個周期內(nèi)的浮子周圍渦量場圖，不難看出，迎著波浪的這一側(cè)會產(chǎn)生較大的渦，因為受到浮子阻礙的流體與波浪傳遞過來的流體之間產(chǎn)生了較大的速度梯度，從而產(chǎn)生了渦。圖中，浮子錐底角度越大，附近渦量分布越紊亂，渦分布的范圍越大，其中90°錐底浮子渦分布范圍很小，僅在浮子下方有一部分，而180°平底浮子流場渦分布范圍最廣。這也就會導(dǎo)致波浪能能量的耗散增加，影響浮子對能量的吸收，進而影響效率?？梢钥闯?，主要產(chǎn)生渦的部位是圓柱轉(zhuǎn)為圓錐處的拐角，這是因為幾何尺寸突然變化，導(dǎo)致流體的速度梯度也變大，容易產(chǎn)生渦。并且錐底角度越大，其尺寸變化幅度也就越大，所產(chǎn)生的大量級的渦也越多，這說明波浪能量通過浮子壁面將能量傳遞給了流體質(zhì)點，變成了流體質(zhì)點的動能。意味著運動過程中耗散的能量增加，這些耗散的能量無法被吸收，而入射波能量固定，這就導(dǎo)致錐底角度小的浮子能量利用率高于角度大的浮子。

圖13 各浮子在t=NT+T/4時的渦量場

圖14 各浮子在t=NT+T/2時的渦量場

圖15 各浮子在t=NT+T3/4時的渦量場

圖16 各浮子在t=(N+1)T時的渦量場

5 結(jié)論

1) 對于平底圓柱浮子，研究了不同直徑吃水比的粘性阻尼大小及水動力性能，圓柱浮子越矮胖，其所受粘性阻尼越小，并且隨著2a/d的增大，浮子俘獲寬度比的峰值區(qū)域也會改變。

2) 對于錐底圓柱浮子，分析了不同錐底角度浮子的水動力性能，當(dāng)錐底角度較小時，更有利于浮子對波浪能量的俘獲。

3) 通過對浮子進行時域流場分析，得到錐底角度大的浮子容易產(chǎn)生較大范圍及量級的渦，影響波浪能量的吸收,與頻域分析結(jié)果一致。