寧德志，石進(jìn)，滕斌，趙海濤

(1.大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧大連116024;2.國(guó)家海洋局第二海洋研究所工程海洋學(xué)重點(diǎn)實(shí)驗(yàn)室，浙江杭州310012)

海洋面積占地球表面積的71%，它不僅為人 類提供航運(yùn)、水產(chǎn)和豐富的礦藏，而且蘊(yùn)含著大量的能量。海洋能主要包括波浪能，潮汐能，海洋溫差能，海洋鹽差能和海流能等，其中波浪能作為一種重要的可再生清潔能源具有分布廣泛、能量密度大、采集結(jié)構(gòu)簡(jiǎn)單等特點(diǎn)。振蕩水柱式波能轉(zhuǎn)換裝置(oscillating water column，OWC)是目前各國(guó)最為重視、公認(rèn)的最有前途、投入研究力量最大、建造使用最多的一種裝置。該裝置通過(guò)氣室將波浪能轉(zhuǎn)換為空氣氣流的能量，再通過(guò)空氣透平將氣流能量轉(zhuǎn)換為電機(jī)轉(zhuǎn)軸的軸功，最后通過(guò)發(fā)電機(jī)將轉(zhuǎn)軸軸功轉(zhuǎn)換為電能。岸式的振蕩水柱轉(zhuǎn)換裝置是振蕩水柱波能轉(zhuǎn)換裝置的一種，它具有結(jié)構(gòu)簡(jiǎn)單，易于安裝和維護(hù)等特點(diǎn)。關(guān)于 OWC的研究有很多。Evans等［1-6］在 OWC的二維和三維頻域勢(shì)流理論分析與數(shù)值計(jì)算中做了很多工作。時(shí)域勢(shì)流理論方面，Koo等［7-8］采用考慮完全非線性自由水面邊界條件的常數(shù)邊界元方法并在氣室內(nèi)引入壓強(qiáng)模型計(jì)算了一種帶底坡的岸式的振蕩水柱波能轉(zhuǎn)換裝置。Zhang等［9-11］應(yīng)用基于N-S方程的粘性流模型對(duì)振蕩水柱式波能轉(zhuǎn)換裝置進(jìn)行了數(shù)值模擬。此外，Sarmento等［12-14］實(shí)驗(yàn)研究了振蕩水柱式波能轉(zhuǎn)換裝置。

本文應(yīng)用時(shí)域高階邊界元方法與壓強(qiáng)模型建立了岸式振蕩水柱波能轉(zhuǎn)換裝置的二維數(shù)值模型，研究了前墻尺寸，氣室寬度及入射波要素對(duì)岸式振蕩水柱波能轉(zhuǎn)換裝置能量轉(zhuǎn)換效率的影響。

1 OWC時(shí)域數(shù)值模型

1.1 控制方程和邊界條件

本文研究的岸式振蕩水柱波能轉(zhuǎn)換裝置如圖1所示，H代表水深，A、B、C、D分別代表氣孔寬度、氣室寬度、前墻壁厚和入水深度。在流體無(wú)粘、無(wú)旋、不可壓縮的假定下，計(jì)算域內(nèi)流體的速度可用速度勢(shì)的梯度表示:

圖1 岸式振蕩水柱波能轉(zhuǎn)換裝置示意圖Fig.1 Schematic of land-based OWC

本模型采用源造波方法造波，控制方程為泊松方程:

式中:q*(xs，z，t)=2vδ(x-xs)為源造波強(qiáng)度，v為流體質(zhì)點(diǎn)的水平速度，本文給定為二階Stokes波速度解析解，波浪產(chǎn)生位置x=xs。

在自由水面上，滿足完全非線性運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)邊界條件:

式中:η為自由水面的鉛垂位移，P為自由水面上的壓強(qiáng)，g為重力加速度。本模型采用混合歐拉-拉格朗日法更新自由水面邊界條件，利用物質(zhì)導(dǎo)數(shù)公式d/dt=?/?t+v·?，并在源造波上游區(qū)域的自由水面布置人工阻尼層防止波浪反射，自由水面邊界條件轉(zhuǎn)化成以下形式:

其中:

式中:X0=(x0，0)為水質(zhì)點(diǎn)初始位置，造波位置為x=0，α為阻尼層的阻尼系數(shù)，βλ為阻尼層的長(zhǎng)度，λ為入射波波長(zhǎng)，β為阻尼層寬度系數(shù)。

在岸式振蕩水柱波能轉(zhuǎn)換裝置表面、水槽側(cè)壁和水底等固定物體邊界上，流體的法向速度為0，即

本模型在時(shí)域內(nèi)求解，在靜水面上初始條件為

1.2 氣室內(nèi)的壓強(qiáng)模型

振蕩水柱式波能轉(zhuǎn)換裝置氣室內(nèi)的氣體壓強(qiáng)與液面的運(yùn)動(dòng)互相影響，存在復(fù)雜的氣液耦合現(xiàn)象，需要引入壓強(qiáng)模型加以處理。本模型假定氣室內(nèi)的壓強(qiáng)與孔口處的氣體流速呈線性關(guān)系［7］:

或二次關(guān)系［7］:

式中:Ud(t)表示氣孔處的氣體流速，Cdm和Ddm分別為線性阻尼系數(shù)和二次阻尼系數(shù)，它們?nèi)Q于氣孔處透平的物理特性，例如威爾斯透平兩端的氣體壓強(qiáng)差與流速就基本滿足線性關(guān)系，式(8)中的絕對(duì)值符號(hào)是為了保證壓強(qiáng)符號(hào)和孔口處氣體流速方向的關(guān)聯(lián)性。在孔口氣體流速與氣室內(nèi)壓強(qiáng)很小的情況下可以假設(shè)空氣不可壓縮，由此可得

1.3 能量轉(zhuǎn)換效率的計(jì)算

入射波浪會(huì)使振蕩水柱波能轉(zhuǎn)換裝置氣室內(nèi)的水體做上下振動(dòng)，從而推動(dòng)氣室內(nèi)的氣體在氣孔兩側(cè)往復(fù)運(yùn)動(dòng)，推動(dòng)透平做功。波浪推動(dòng)氣體做功的速率可由下式求得

仿照波能流的推導(dǎo)，在一個(gè)周期內(nèi)做平均可得

式中:z(t)表示氣室內(nèi)平均水面的速度。假定氣室內(nèi)的空氣不可壓縮且孔口處的氣體流速呈正弦變化:

將式(8)、(12)代入式(11)可以推得

氣室內(nèi)的能量轉(zhuǎn)換效率κ可由下式求得

式中:分母為入射波波能流，Ai為入射波波幅，Cg為群速度。

1.4 數(shù)值求解

在流體計(jì)算域內(nèi)對(duì)速度勢(shì)應(yīng)用第二格林公式，可得到如下邊界積分方程［15］:

式中:p=(x0，z0)為源點(diǎn)，q=(x，z)為場(chǎng)點(diǎn)，C(p)為固角系數(shù)，Γ為包括自由水面邊界和固體邊界的流體計(jì)算邊界，Ω代表整個(gè)流體計(jì)算域。為了消除水底的積分，減小計(jì)算量，選取考慮水底鏡像的格林函數(shù):

式中:r1為p和q兩點(diǎn)距離，r2為p和q關(guān)于水底的鏡像之間距離。

本文用三節(jié)點(diǎn)的二次單元將計(jì)算域邊界離散成一些曲線單元，每個(gè)單元可以通過(guò)數(shù)學(xué)變換變換成參數(shù)坐標(biāo)下的等參元，在等參元內(nèi)引入形函數(shù)可以插值得到單元內(nèi)任一點(diǎn)的幾何坐標(biāo)和速度勢(shì)。將源點(diǎn)取在各個(gè)節(jié)點(diǎn)上可以將離散后的積分方程化為線性方程組的形式，從而求解未知量。每一步的計(jì)算都認(rèn)為當(dāng)前時(shí)刻自由水面上的速度勢(shì)和物面上的速度勢(shì)法向?qū)?shù)是已知的，通過(guò)積分方程計(jì)算當(dāng)前時(shí)刻自由水面上的速度勢(shì)法向?qū)?shù)和物面上的速度勢(shì)。在氣室內(nèi)應(yīng)用壓強(qiáng)模型后，下一時(shí)刻的自由水面上水質(zhì)點(diǎn)的位置和速度勢(shì)可通過(guò)四階Runga-Kutta法計(jì)算，然后重新劃分網(wǎng)格，繼續(xù)應(yīng)用積分方程計(jì)算自由水面上的速度勢(shì)法向?qū)?shù)和物面上的速度勢(shì)。這樣計(jì)算周而復(fù)始，直到計(jì)算結(jié)束［16］。

2 OWC數(shù)值模型驗(yàn)證

氣室內(nèi)的壓強(qiáng)模型已經(jīng)成功被用于氣動(dòng)浮式防波堤的研究中，為了驗(yàn)證本文壓強(qiáng)模型的準(zhǔn)確性，對(duì)比計(jì)算了文獻(xiàn)［7］研究過(guò)的一種氣動(dòng)防波堤的相應(yīng)工況。計(jì)算水深30 m，氣室寬度30 m，前、后墻壁厚10 m，入水深度10 m，氣孔寬度1 m。入射波周期9～12 s，入射波波幅0.5 m。本研究中計(jì)算域長(zhǎng)度取10倍波長(zhǎng)距離，氣室箱體布置在計(jì)算域中間，水槽兩端均布置長(zhǎng)度為1.5λ的阻尼層。通過(guò)數(shù)值收斂性試驗(yàn)，自由水面上每個(gè)波長(zhǎng)長(zhǎng)度布置20個(gè)網(wǎng)格，沿水深方向布置10個(gè)網(wǎng)格，時(shí)間步長(zhǎng)Δt=T/60。圖2(a)、(b)分別給出了入射波周期為12 s和10 s時(shí)，透射波高增大率γ隨線性阻尼系數(shù)Cdm與二次阻尼系數(shù)Ddm變化的對(duì)比曲線。從圖中可以看出，本文結(jié)果與文獻(xiàn)［7］的數(shù)值結(jié)果吻合很好。

圖2 不同壓強(qiáng)模型的對(duì)比Fig.2 Comparison of the pressure models

圖3是線性壓強(qiáng)模型下，Cdm=150，入射波周期T=10 s時(shí)防波堤氣室內(nèi)壓強(qiáng)的時(shí)間歷程曲線?？梢钥吹綒馐覂?nèi)壓強(qiáng)波動(dòng)周期與入射波周期相同。本模型氣室內(nèi)的壓強(qiáng)在一定時(shí)間后趨于穩(wěn)定，說(shuō)明本模型可以消除二次反射的影響進(jìn)行長(zhǎng)時(shí)間模擬。

圖3 氣動(dòng)浮式防波堤氣室內(nèi)壓強(qiáng)時(shí)間歷程，Cdm=150，T=10 sFig.3 Time history of the pressure in the air chamber of pneumatic floating breakwater，Cdm=150，T=10 s

為了進(jìn)一步驗(yàn)證數(shù)值模型的準(zhǔn)確性，將本文結(jié)果與文獻(xiàn)［3，9，14］解析解、數(shù)值解和實(shí)驗(yàn)值進(jìn)行了對(duì)比。振蕩水柱式波能轉(zhuǎn)換裝置如圖1，具體參數(shù)見(jiàn)表1。氣孔寬度A=0.005 m，計(jì)算水深H=0.92 m，入射波幅0.04 m。計(jì)算域長(zhǎng)度取5倍波長(zhǎng)距離，通過(guò)開(kāi)展數(shù)值收斂性驗(yàn)證，單位波長(zhǎng)長(zhǎng)度內(nèi)劃分30個(gè)網(wǎng)格，時(shí)間步長(zhǎng)Δt=T/80。本文采用線性壓強(qiáng)模型，經(jīng)過(guò)試算，阻尼系數(shù)Cdm選定為 3.5。

表1 岸式振蕩水柱波能轉(zhuǎn)換裝置結(jié)構(gòu)尺寸Table 1 Geometrical parameters of the land-based OWC

圖4給出了3種前墻尺寸的波能轉(zhuǎn)換效率κ隨無(wú)量綱入射波周期T*變化的對(duì)比結(jié)果，T*=從圖中可以看出在高頻區(qū)(T*＜6)，3種不同尺寸的模型計(jì)算中，本文結(jié)果和文獻(xiàn)［9，14］吻合的都較好。在低頻區(qū)(T*＞7)，本文結(jié)果與文獻(xiàn)［9］的數(shù)值結(jié)果吻合較好，但與文獻(xiàn)［14］的實(shí)驗(yàn)結(jié)果相比數(shù)值偏小，可能的原因是在實(shí)驗(yàn)中的造波板與結(jié)構(gòu)距離為37.5 m，在入射波波長(zhǎng)較長(zhǎng)的情況下二次反射很快會(huì)影響到結(jié)構(gòu)處的測(cè)量結(jié)果，造成能量轉(zhuǎn)換效率被高估。本文結(jié)果與文獻(xiàn)［3］解析解的趨勢(shì)和共振頻率基本一致，但數(shù)值上有一定差別，這是由于解析解對(duì)模型做了很多理想化假定，波能的最大轉(zhuǎn)化效率可以達(dá)到100%，這在實(shí)際中是不可能達(dá)到的。總體來(lái)說(shuō)，相對(duì)于文獻(xiàn)［9］的兩相流模型，本模型無(wú)論是可計(jì)算的頻率寬度還是與實(shí)驗(yàn)結(jié)果的對(duì)比都較好，具有較高的計(jì)算精度。

圖4 不同前墻厚度下，能量轉(zhuǎn)換效率隨無(wú)量綱入射周期T*的變化Fig.4 Energy conversion efficiency versus T*with different thickness of the front wall

3 結(jié)構(gòu)尺寸對(duì)能量轉(zhuǎn)換效率的影響

在水深一定的情況下，岸式振蕩水柱波能轉(zhuǎn)換裝置的前墻入水深度、前墻厚度和氣室寬度是影響其能量轉(zhuǎn)換效率的主要因素。本文計(jì)算水深統(tǒng)一為0.92 m，入射波幅0.04 m。圖5給出了在氣室寬0.64 m，前墻厚度0.04 m時(shí)，結(jié)構(gòu)在不同前墻入水深度下的能量轉(zhuǎn)換效率對(duì)比曲線。從圖中可以看到，整個(gè)結(jié)構(gòu)存在一個(gè)明顯的共振頻率，在此頻率上的能量轉(zhuǎn)換效率最高，可以達(dá)到70%以上。隨著入射波頻率遠(yuǎn)離共振頻率，振蕩水柱波能轉(zhuǎn)換裝置的能量轉(zhuǎn)換效率逐漸變小。在當(dāng)前計(jì)算范圍內(nèi)，前墻入水深度的增大會(huì)導(dǎo)致結(jié)構(gòu)的共振頻率區(qū)間向低頻區(qū)域移動(dòng)，同時(shí)高頻區(qū)(T*＜6)的能量轉(zhuǎn)換效率變小，共振頻率上的能量轉(zhuǎn)換效率增大。前墻入水深度的變化對(duì)低頻區(qū)(T*＞7)的能量轉(zhuǎn)換效率影響不大，這是因?yàn)樵诘皖l區(qū)前墻入水深度相對(duì)于波長(zhǎng)很小。

圖6給出了氣室寬度0.64 m，前墻入水深度0.15 m情況下，3種不同前墻厚度下的能量轉(zhuǎn)換效率曲線。同樣，在低頻區(qū)(T*＞7)，前墻厚度的變化對(duì)能量轉(zhuǎn)換效率的影響很小。隨著前墻厚度的增大結(jié)構(gòu)的共振區(qū)間向低頻區(qū)的方向移動(dòng)，在高頻區(qū)和共振頻率上，前墻厚度的增大會(huì)導(dǎo)致能量轉(zhuǎn)換效率的減小。

圖7給出了在前墻入水深度0.15 m，前墻厚度0.04 m，不同氣室寬度下的能量轉(zhuǎn)換效率對(duì)比曲線。改變氣室的寬度會(huì)同時(shí)改變氣室在低頻區(qū)與高頻區(qū)的能量轉(zhuǎn)換效率，隨著氣室寬度的增大結(jié)構(gòu)的共振頻率向低頻方向移動(dòng)，且能量轉(zhuǎn)換效率最大值會(huì)隨氣室寬度的增大而變小。氣室寬度的增大會(huì)改善結(jié)構(gòu)低頻區(qū)的能量轉(zhuǎn)換效率，降低結(jié)構(gòu)在高頻區(qū)的能量轉(zhuǎn)換效率。

圖5 不同前墻入水深度下，能量轉(zhuǎn)換效率隨T*的變化Fig.5 Energy conversion efficiency versus T*with different immersion depth of the front wall

圖6 不同前墻厚度下，能量轉(zhuǎn)換效率隨T*的變化Fig.6 Energy conversion efficiency versus T*with different thickness of front wall

圖7 不同氣室寬度下，能量轉(zhuǎn)換效率隨T*的變化Fig.7 Energy conversion efficiency versus T*with different width of the chamber

4 結(jié)論

本研究利用時(shí)域高階邊界元理論與氣體壓強(qiáng)模型建立了岸式振蕩水柱波能轉(zhuǎn)換裝置的完全非線性數(shù)值模型。從研究中可以發(fā)現(xiàn)岸式振蕩水柱波能轉(zhuǎn)換裝置具有一個(gè)明顯的共振區(qū)間，適當(dāng)減小前墻入水深度、厚度與氣室寬度都可以增大共振頻率、提高高頻區(qū)的能量轉(zhuǎn)換效率;前墻尺寸的變化對(duì)低頻區(qū)能量轉(zhuǎn)換效率的影響十分微小;氣室寬度的增加會(huì)提高低頻區(qū)的能量轉(zhuǎn)換效率，降低高頻區(qū)的能量轉(zhuǎn)換效率。和前人的勢(shì)流模型相比，本模型采用高階單元與源造波技術(shù)，精度更高，所需單元更少且消波問(wèn)題很容易處理。和粘性流模型相比，本模型的計(jì)算效率更高，在合理的阻尼系數(shù)下同樣可以與實(shí)驗(yàn)結(jié)果吻合良好。由于無(wú)法計(jì)算由于流動(dòng)分離、渦旋等引起的粘性耗散并且尚未考慮空氣的可壓縮性，本模型存在一定的局限性，下一步工作需要建立更加貼近實(shí)際的壓強(qiáng)模型并考慮粘性耗散的影響，同時(shí)對(duì)模型進(jìn)行三維拓展，以達(dá)到更好的模擬效果。

［1］EVANS D.The oscillating water column wave-energy device［J］.IMA Journal of Applied Mathematics，1978，22(4):423-433.

［2］游亞戈.復(fù)雜地形下岸式波能裝置的水動(dòng)力學(xué)計(jì)算［J］.海洋技術(shù)，1993，12(1):20.YOU Yage.Hydrodynamic calculation of shore-mounted wave-power device at complicated topography［J］.Ocean Tecnology，1993，12(1):20.

［3］EVANS D，PORTER R.Hydrodynamic characteristics of an oscillating water column device［J］.Applied Ocean Research，1995，17(3):155-164.

［4］DELAUR Y，LEWIS A.3D hydrodynamic modelling of fixed oscillating water column wave power plant by a boundary element methods［J］.Ocean Engineering，2003，30(3):309-330.

［5］WANG D，KATORY M，LI Y.Analytical and experimental investigation on the hydrodynamic performance of onshore wave-power devices［J］.Ocean Engineering，2002，29(8):871-885.

［6］HONG D，HONG S Y，HONG S W.Numerical study of the motions and drift force of a floating OWC device［J］.Ocean Engineering，2004，31(2):139-164.

［7］KOO W，KIM M，LEE D.Nonlinear time-domain simulation of pneumatic floating breakwater［J］.International Journal of Offshore and Polar Engineering，2006，16(1):25-32.

［8］KOO W，KIM M-H.Nonlinear time-domain simulation of a land-based oscillating water column［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2010，136(5):276-285.

［9］ZHANG Y，ZOU Q P，GREAVES D.Air-water two-phase flow modelling of hydrodynamic performance of an oscillating water column device［J］.Renewable Energy，2011，41:159-170.

［10］劉臻.岸式振蕩水柱波能發(fā)電裝置的試驗(yàn)及數(shù)值模擬研究［D］.青島:中國(guó)海洋大學(xué)，2008:70-107.LIU Zhen.Experimental and numerical investigation of oscillating water column wave energy converter［D］.Qingdao:Ocean University of China，2008:70-107.

［11］紀(jì)君娜，劉臻，紀(jì)立強(qiáng).振蕩水柱波能發(fā)電裝置氣室的三維數(shù)值模擬研究［J］.海岸工程，2011，30(2):7-13.JI Junna，LIU Zhen，JI Liqiang.3D numerical simulation for oscillating water column chamber in wave-power converter［J］.Coastal Engineering，2011，30(2):7-13.

［12］SARMENTO A，F(xiàn)ALCAO A F O.Wave generation by an oscillating surface-pressure and its application in wave-energy extraction［J］.Journal of Fluid Mechanics，1985，150:467-485.

［13］TSENG R S，WU R H，HUANG C C.Model study of a shoreline wave-power system［J］.Ocean Engineering，2000，27(8):801-821.

［14］MORRIS-THOMAS M T，IRVIN R J，THIAGARAJAN K P.An investigation into the hydrodynamic efficiency of an oscillating water column［J］.Journal of Offshore Mechanics and Arctic Engineering，2007，129(4):273-278.

［15］NING D，TENG B，EATOCK T R，et al.Numerical simulation of non-linear regular and focused waves in an infinite water-depth［J］.Ocean Engineering，2008，35(8):887-899.

［16］NING D，TENG B.Numerical simulation of fully nonlinear irregular wave tank in three dimension［J］.International Journal for Numerical Methods in Fluids，2007，53(12):1847-1862.