劉葳興，劉 磊，陳雨龍，王 鑫，張之陽(yáng)

（江蘇海洋大學(xué)海洋工程學(xué)院，江蘇 連云港 222005）

潮汐能對(duì)比其他可再生能源具有再生、可預(yù)測(cè)、環(huán)保等優(yōu)勢(shì)[1]，從潮汐能中提取能量是當(dāng)下研究熱點(diǎn)，而潮汐能水輪機(jī)是一種從潮汐能中提取能量轉(zhuǎn)換為其他能量的裝置。模型試驗(yàn)、數(shù)值仿真可用于水輪機(jī)的尾流場(chǎng)研究[2-3]。水輪機(jī)的數(shù)值仿真方法常分為全尺寸仿真、縮放尺寸仿真和致動(dòng)盤(pán)方法。然而，從節(jié)約計(jì)算時(shí)間考慮，可用使致動(dòng)盤(pán)方法捕捉水輪機(jī)工作時(shí)尾流。用致動(dòng)盤(pán)方法替代真實(shí)葉片進(jìn)行計(jì)算時(shí)，常用體積力作為附加源項(xiàng)，并規(guī)定體積力分布形式及其大小[4]。致動(dòng)盤(pán)方法在捕捉葉片幾何形狀變化對(duì)流場(chǎng)造成的影響時(shí)有所欠缺。捕捉翼型對(duì)流場(chǎng)帶來(lái)的影響可使用葉片元素動(dòng)量理論（blade element momentum theory，BEM）方法[5-6]。張之陽(yáng)等[7]采用BEM 分析水平軸潮流能渦輪機(jī)的水動(dòng)力特性，其方法節(jié)省大量計(jì)算時(shí)間，但在尾流細(xì)節(jié)處研究不足。而應(yīng)用致動(dòng)盤(pán)模型（actua‐tor disc model），可在計(jì)算時(shí)間上、節(jié)省計(jì)算資源和捕捉尾流細(xì)節(jié)特征上相平衡。張玉全等[8]基于曼徹斯特大學(xué)循環(huán)水槽的實(shí)驗(yàn)數(shù)據(jù)，利用致動(dòng)盤(pán)理論求解Navier-Stokes 方程獲得尾流場(chǎng)數(shù)據(jù)，并將兩者對(duì)比，發(fā)現(xiàn)致動(dòng)盤(pán)方法更易獲得尾流場(chǎng)情況。Myers等[9]用多孔圓盤(pán)代替潮流能水輪機(jī)葉片，對(duì)其性能開(kāi)展研究，發(fā)現(xiàn)水輪機(jī)的推力、機(jī)架離拖曳池底面距離與水池壁面粗糙度對(duì)尾流都有不同影響。Sun等[10]使用Fluent軟件建立多孔圓盤(pán)的二維和三維模型，發(fā)現(xiàn)水輪機(jī)運(yùn)行中流場(chǎng)后方的液面高度有所下降與尾流的速度虧損有關(guān)。Edmunds等[11]使用非線性致動(dòng)盤(pán)模型并對(duì)傳統(tǒng)源項(xiàng)進(jìn)行改進(jìn)，分析葉片尖端損失及尾流場(chǎng)。致動(dòng)盤(pán)替代實(shí)際葉片進(jìn)行仿真，常用于風(fēng)機(jī)尾流場(chǎng)分析，因風(fēng)機(jī)與水平軸水輪機(jī)在結(jié)構(gòu)上具有一定相似性，故可將致動(dòng)盤(pán)模型應(yīng)用于水輪機(jī)尾流場(chǎng)分析。本研究采用體積力附加源項(xiàng)方式對(duì)潮汐能水輪進(jìn)行數(shù)值仿真，并與愛(ài)丁堡大學(xué)中的拖曳池試驗(yàn)作對(duì)比，評(píng)估致動(dòng)盤(pán)模型在水輪機(jī)上的應(yīng)用。

1 實(shí)驗(yàn)裝置與方法

1.1 試驗(yàn)拖曳池

本實(shí)驗(yàn)數(shù)據(jù)參考愛(ài)丁堡大學(xué)Flowave 拖曳池中進(jìn)行的水輪機(jī)尾流場(chǎng)實(shí)驗(yàn)[12]。實(shí)驗(yàn)拖曳池為直徑25 m 的圓形結(jié)構(gòu)，工作水深H=2 m。拖曳池中，沿整個(gè)水池壁面有28個(gè)葉輪單元，每個(gè)單元可控制水流流速而形成一個(gè)循環(huán)流動(dòng)系統(tǒng)。通過(guò)設(shè)置葉輪轉(zhuǎn)向在水池中央試驗(yàn)區(qū)的任一方向產(chǎn)生水流流動(dòng)，可實(shí)現(xiàn)最大水流流速為1.6 m/s。水池上方有XY 行車(chē)方便安裝水輪機(jī)。由于實(shí)驗(yàn)中沒(méi)有設(shè)施可控制湍流強(qiáng)度，因此，本實(shí)驗(yàn)湍流強(qiáng)度與愛(ài)丁堡大學(xué)進(jìn)行水輪機(jī)實(shí)驗(yàn)所選擇的湍流強(qiáng)度相同（該湍流強(qiáng)度為歐洲海洋平均流速的7%[13]）。

1.2 渦輪機(jī)和儀器設(shè)備

固定在拖曳池底部的水平軸水輪機(jī)，葉片直徑D=1.2 m，輪轂中心距水池底h=1 m。其翼型為NACA63812、葉片的幾何形狀可參考文獻(xiàn)[13]。水輪機(jī)上有安裝測(cè)量葉片彎矩R，轉(zhuǎn)子扭矩Q和推力T的傳感器。水輪機(jī)中的伺服電機(jī)可使用編碼器來(lái)控制轉(zhuǎn)子的轉(zhuǎn)動(dòng)角度θ，進(jìn)而使葉片在不同尖端速度比（tip-speed ratio，TSR）范圍內(nèi)運(yùn)行。

水輪機(jī)流場(chǎng)布置見(jiàn)圖1（a），拖曳池水流方向沿x軸，水池縱向?yàn)閥軸，深度方向?yàn)閦軸，轉(zhuǎn)子中心位置記（x,y,z）=（0,0,0）。用Vectrino Profiler聲學(xué)多普勒測(cè)速儀（ADV）對(duì)流量采樣，以轉(zhuǎn)子中心為水平面（z=1 m 截面），采樣點(diǎn)在x軸、x=0 m、x=2.4 m 上分布，見(jiàn)圖1（b）。采樣頻率為100 Hz，采樣時(shí)間256 s。測(cè)量的流量速度分別為水平速度、縱向速度和垂直速度。實(shí)驗(yàn)所用水流速度為0.8 m/s。

圖1 水輪機(jī)布置示意Fig.1 Schematic showing the turbine arrangement

水輪機(jī)實(shí)驗(yàn)時(shí)還需驗(yàn)證其運(yùn)行過(guò)程中是否產(chǎn)生阻塞效應(yīng)。水輪機(jī)阻塞率為葉片旋轉(zhuǎn)一周面積S1比拖曳池面積S2的值，S1/S2=0.002 3。垂直阻塞率為葉片直徑D比水池水深h的值，D/h=0.6。水平阻塞率為葉片直徑D比水池直徑w的值，D/w=0.048?？梢?jiàn)，水輪機(jī)阻塞率、垂直阻塞率、水平阻塞率均很小，表明在潮汐能水輪機(jī)測(cè)試期間不會(huì)出現(xiàn)阻塞效應(yīng)進(jìn)而影響流場(chǎng)。

2 致動(dòng)盤(pán)模型

致動(dòng)盤(pán)方法基于一維動(dòng)量理論，將水輪機(jī)葉輪視為一個(gè)可穿透的等效面積圓盤(pán)區(qū)域代替[14]。同時(shí)，將水輪機(jī)致動(dòng)盤(pán)周?chē)鲌?chǎng)視為穩(wěn)態(tài)、一維、不可壓流場(chǎng)，并以體積力的形式實(shí)現(xiàn)致動(dòng)盤(pán)對(duì)流場(chǎng)的作用[15]。

由于致動(dòng)盤(pán)模型常用于風(fēng)機(jī)尾流場(chǎng)研究，而在水輪機(jī)上應(yīng)用較少，因此，本研究基于文獻(xiàn)[16]，將OpenFOAM 的自帶“simpleFoam”求解器進(jìn)行修改自定義成新的致動(dòng)盤(pán)求解器，并將fan 邊界條件在OpenFOAM 中應(yīng)用。假設(shè)在不可逆流體下，方程組可寫(xiě)成以下形式：

其中，xi代表笛卡爾坐標(biāo)分量（x,y,z），Ui是笛卡爾平均速度分量，fi包括圓盤(pán)轉(zhuǎn)子特征的附加源，即葉片對(duì)流體作用，本研究所用體積力作為附加源項(xiàng)。雷諾應(yīng)力為，需應(yīng)用湍流模型來(lái)建模得到控制方程。

將實(shí)驗(yàn)得到的推力T以體積力源項(xiàng)形式均勻分布于致動(dòng)盤(pán)模型上。體積力與來(lái)流速度和推力系數(shù)有關(guān)。假設(shè)來(lái)流速度為U，對(duì)應(yīng)功率曲線的空氣密度和推力系數(shù)分別為ρ,CT，則水輪機(jī)葉輪單位面積的軸向力為

Δx為致動(dòng)盤(pán)厚度，V為致動(dòng)盤(pán)體積。其動(dòng)量方程中，致動(dòng)盤(pán)體積力源項(xiàng)表示為[17]

本研究采用RNGk-ε湍流模型，k代表湍流動(dòng)能，ε代表這個(gè)能量的耗散率。

3 數(shù)值模擬

3.1 網(wǎng)格規(guī)劃

仿真建模時(shí)對(duì)拖曳池和水輪機(jī)按1∶1 建模，以保證數(shù)值仿真的準(zhǔn)確性。用OpenFOAM 6.0版本軟件中“blockMesh”和“snappyHexMesh”功能劃分網(wǎng)格?！癰lockMesh”工具用來(lái)生成一個(gè)初始?jí)K（網(wǎng)格域），并將其細(xì)分為不同的尺寸區(qū)域。塊初始尺寸的底部面積為30 m2，高4 m 的圓柱體。接著“simpleGrading”參數(shù)設(shè)置為1。生成塊狀網(wǎng)格后，使用“snappyHexMesh”工具迭代細(xì)化“blockMesh”，將產(chǎn)生的拆分的“hex”網(wǎng)格變形到幾何體上，從而近似地符合幾何體。

將尾流區(qū)域規(guī)劃為一個(gè)半徑為0.7 m 的圓柱體，尾流長(zhǎng)度以轉(zhuǎn)子中心起到尾流下游9 m 處結(jié)束。對(duì)流體計(jì)算域不同區(qū)域的網(wǎng)格進(jìn)行加密（圖2），水輪機(jī)周?chē)吔鐚泳W(wǎng)格密度，細(xì)化水平最高；機(jī)架上方網(wǎng)格密度次之；為節(jié)省計(jì)算資源，流體計(jì)算外域網(wǎng)格密度最低。計(jì)算域網(wǎng)格規(guī)劃完成后進(jìn)行網(wǎng)格敏感性分析，網(wǎng)格數(shù)量分4 個(gè)等級(jí)（表1）。對(duì)上述4 種數(shù)目的網(wǎng)格進(jìn)行計(jì)算，計(jì)算結(jié)果以y=0 m，z=1 m截面速度分析，計(jì)算結(jié)果見(jiàn)圖3。由圖3 可見(jiàn)，當(dāng)選用網(wǎng)格數(shù)量等級(jí)3后，流場(chǎng)中速度差異變化較小，故為節(jié)省計(jì)算時(shí)間后續(xù)仿真計(jì)算選擇網(wǎng)格數(shù)量等級(jí)3進(jìn)行求解。

圖2 渦輪機(jī)網(wǎng)格細(xì)化示意Fig.2 Turbine mesh refinement diagram

表1 網(wǎng)格數(shù)量Table 1 Total number of meshes

圖3 網(wǎng)格敏感度Fig.3 Mesh sensitivity

3.2 初始條件和邊界條件

水池壁面各部分按不同位置設(shè)置為質(zhì)量流入、質(zhì)量流出或無(wú)流流動(dòng)。將圖4 拖曳池中的28 個(gè)葉輪單元配置為入口、出口及壁面，質(zhì)量流量大小根據(jù)表2 設(shè)定。表2 中質(zhì)量流量（X軸流向）正值為流入，負(fù)值為流出。流入之和與流出之和為相同數(shù)值，則保持質(zhì)量守恒。水池壁面速度設(shè)置為0 m/s，壁面函數(shù)用于k，ε。流體域頂部設(shè)置為無(wú)滑移壁面。除速度外，其余初始條件皆設(shè)置為邊界條件，初始速度設(shè)置為0 m/s，動(dòng)力粘度ν設(shè)定為1.666 7 ×10-6m2/s?？刂品匠糖蠼膺x擇以壓力為基礎(chǔ)的求解器，采用Simple 算法，空間離散采用二階迎風(fēng)格式，收斂標(biāo)準(zhǔn)為1 × 10-4。

圖4 拖曳池質(zhì)量流量配置示意Fig.4 Mass flow configuration diagram for towing pool

表2 質(zhì)量流量設(shè)置Table 2 Mass flow rate setting

4 仿真結(jié)果與分析

本研究仿真結(jié)果分為兩部分，一是拖曳池?cái)?shù)值仿真，即在無(wú)水輪機(jī)情況下水流在拖曳池中自由發(fā)展；二是在致動(dòng)盤(pán)模型下進(jìn)行仿真。將兩者仿真結(jié)果對(duì)比，分析水輪機(jī)尾流虧損情況，并為后續(xù)研究做基礎(chǔ)。

4.1 自由發(fā)展區(qū)仿真

以水輪機(jī)轉(zhuǎn)子中心z=1 為截面拖曳池?cái)?shù)值仿真見(jiàn)圖5。當(dāng)水流速度以0.8 m/s 從壁面入口加速流向拖曳池中心時(shí)，水流會(huì)與壁面發(fā)生相互作用，引起湍流強(qiáng)度的加劇，從而導(dǎo)致速度虧損。水流在拖曳池中心充分發(fā)展時(shí)，流動(dòng)速度更加均勻。水流經(jīng)下游壁面流出時(shí)，部分水流相互作用時(shí)會(huì)產(chǎn)生出口回流，導(dǎo)致速度在壁面附近發(fā)散，又因水池上下壁面不是出入口，回流在此加劇湍流強(qiáng)度。

圖5 拖曳池仿真Fig.5 Simulation diagram of fluid flow rate in the towing pool

以圖1（b）采樣點(diǎn)對(duì)拖曳池不同位置的流速和湍流數(shù)據(jù)進(jìn)行采樣評(píng)估（圖6）。由圖6（a）可見(jiàn)，上游水流速度沿x軸逐漸增大，到x=0 附近水流增速變緩，x=0 之后速度開(kāi)始下降。這是由于上下游水流與壁面作用引起速度虧損，以及拖曳池上下湍流強(qiáng)度大引起水流流速在中心處增大和水流經(jīng)充分發(fā)展在拖曳池中心處流速均勻。圖6（a）中，實(shí)驗(yàn)與仿真流速誤差在靠近x=0位置時(shí)減少，并隨著尾流長(zhǎng)度的增加而大，湍流強(qiáng)度兩者誤差較小。實(shí)驗(yàn)和仿真計(jì)算的水流速度和湍流強(qiáng)度之間的平均誤差分別為0.3%和1.2%。

圖6（b，c）顯示以x軸向上的水流發(fā)展和擴(kuò)散，在仿真模擬結(jié)果中可見(jiàn)速度和湍流強(qiáng)度曲線的對(duì)稱(chēng)性，然而在實(shí)驗(yàn)中存在不對(duì)稱(chēng)性。這可能是在實(shí)驗(yàn)過(guò)程中很難精確控制拖曳池中各個(gè)單元所需的水流流量，從而引起整個(gè)區(qū)域的平均流量和湍流強(qiáng)度變化。在x=0 m 和x=2.4 m處，實(shí)驗(yàn)和仿真數(shù)值的水流速度的差值分別為6.7%和8.2%。兩截面還顯示出外側(cè)水流速度均大于中心水流速度，其原因是拖曳池外側(cè)葉輪控制的質(zhì)量流入大于中心處的質(zhì)量流入。

圖6 不同截面速度與湍流強(qiáng)度Fig.6 Different cross-sectional velocities and turbulence intensities

4.2 致動(dòng)盤(pán)模型仿真

控制水輪機(jī)葉片轉(zhuǎn)速，進(jìn)而改變尖端速度比（TSR），實(shí)驗(yàn)中TSR 設(shè)置為4～7。由于推力T以體積力源項(xiàng)形式均勻分布于致動(dòng)盤(pán)模型上，因此需要得到推力系數(shù)CT。由前文實(shí)驗(yàn)，得到推力T，帶入公式（5），便可得到推力系數(shù)CT。不同的TSR 可以評(píng)估功率系數(shù)CP，進(jìn)而得到水輪機(jī)最大功率。TSR 是由葉片旋轉(zhuǎn)速度比來(lái)流速度計(jì)算得來(lái)，見(jiàn)公式（6）。CP由水輪機(jī)實(shí)際功率比上理論功率得來(lái)，見(jiàn)公式（7）。由圖7 可見(jiàn)，TSR=6 時(shí)，功率系數(shù)CP最大，且只比TSR=5.5 時(shí)的CP略大，又因TSR=5.5 時(shí)實(shí)驗(yàn)和仿真的推力系數(shù)CT高度吻合，故綜合考慮后續(xù)仿真實(shí)驗(yàn)采取TSR=5.5進(jìn)行。

圖7 不同TSR下CP與CT值的比較Fig.7 CP and CT results under a range of TSRs

其中，T為推力，ρ流體密度,U為來(lái)流速度。

其中，ω為旋轉(zhuǎn)角速度。

其中，P為水輪機(jī)實(shí)際功率，A為致動(dòng)盤(pán)面積。

圖8（a）為z=1 m 截面水輪機(jī)速度云圖，左圖為速度細(xì)節(jié)，右圖為流場(chǎng)區(qū)整體速度。圖8（b）為湍流強(qiáng)度云圖，左圖為水輪機(jī)尾流區(qū)湍流強(qiáng)度細(xì)節(jié)，右圖為流場(chǎng)區(qū)整體湍流強(qiáng)度。圖10為y=0 m，z=1 m截面水輪機(jī)速度和速度虧損曲線。由圖9 可見(jiàn)，速度和湍流數(shù)據(jù)的實(shí)驗(yàn)值與數(shù)值模擬高度吻合，兩者誤差分別為3%到5%。結(jié)合圖8（a）和圖9 可看出，水流速度由上游到致動(dòng)盤(pán)前1.5 m（x=-1.5 m）處增加，接著由于轉(zhuǎn)動(dòng)盤(pán)及機(jī)架對(duì)水流有阻礙作用引起湍流而使水流速度下降。在致動(dòng)盤(pán)上存在不連續(xù)的壓強(qiáng)下降，致使尾流向兩側(cè)擴(kuò)張。水流在經(jīng)過(guò)旋轉(zhuǎn)的轉(zhuǎn)盤(pán)后，由于轉(zhuǎn)盤(pán)吸收來(lái)流的部分能量并且來(lái)流與機(jī)架相互作用，使尾流動(dòng)能小于來(lái)流動(dòng)能，該區(qū)域內(nèi)產(chǎn)生較大的湍流，同時(shí)降低流速并且尾流恢復(fù)嚴(yán)重不足。圖8（a）顯示，下游尾流的速度云圖中渦流呈不對(duì)稱(chēng)狀態(tài)，其原因可能是水輪機(jī)建模細(xì)節(jié)處網(wǎng)格劃分不一致。圖9 流量虧損顯示，上游流量到制動(dòng)盤(pán)前，水輪機(jī)與自由發(fā)展區(qū)要之間約有-0.018～0.13 m/s的流量虧損。

圖8 水輪機(jī)仿真Fig.8 Tidal turbine simulation contour

圖9 y=0 m,z=1 m截面速度Fig.9 y=0 m,z=1 m section velocity

觀察尾流流場(chǎng)不同剖面圖的速度變化是衡量尾流水動(dòng)力特性的重要依據(jù)。圖10 為x=-0.5、0、0.5、2.4、4.0 m，z=1 m 處截面圖。其中x=2.4 m 處的實(shí)驗(yàn)數(shù)據(jù)和仿真結(jié)果高度一致，表明致動(dòng)盤(pán)方法在捕捉尾流特性效果較好。從圖10可知，水輪機(jī)從水流中提取動(dòng)能，通過(guò)水輪機(jī)的部分流體會(huì)有減速效應(yīng)，并產(chǎn)生一個(gè)對(duì)稱(chēng)“W 形尾流”。當(dāng)水流在尾流中進(jìn)一步擴(kuò)散時(shí)，“W 形尾流”剖面消失。隨著尾流的擴(kuò)大，尾流剖面寬度將大于水輪機(jī)直徑，并且由于湍流和機(jī)架作用尾流中心處的速度逐漸下降，而兩側(cè)尾流速度逐漸恢復(fù)。這說(shuō)明尾流場(chǎng)兩側(cè)與周?chē)h(huán)境流場(chǎng)的動(dòng)量交換促使兩側(cè)速度恢復(fù)，而中心區(qū)域隨著與湍流區(qū)越接近使其耗散的能量也就越大。

圖10 TSR=5.5下不同截面速度Fig.10 TSR=5.5,different section velocities

5 結(jié)論

本研究應(yīng)用實(shí)驗(yàn)和數(shù)值仿真的方法對(duì)比分析潮汐能水輪機(jī)尾流場(chǎng)特性。為保證仿真結(jié)果準(zhǔn)確性，建模時(shí)將水輪機(jī)機(jī)架結(jié)構(gòu)一并考慮進(jìn)去。在實(shí)驗(yàn)設(shè)計(jì)過(guò)程中，通過(guò)拖曳池中有無(wú)水輪機(jī)對(duì)其流場(chǎng)影響分析。得到以下結(jié)論：

1）拖曳池?cái)?shù)值仿真結(jié)果得到的速度與湍流強(qiáng)度與實(shí)驗(yàn)數(shù)據(jù)相吻合，但在縱向截面（z=1 m，不同y截面）上仿真模擬和實(shí)驗(yàn)之間有輕微的差異。這是因?yàn)樵趯?shí)驗(yàn)過(guò)程中，葉輪機(jī)械存在摩擦、阻力等因素，無(wú)法精確控制水流的輸出與仿真過(guò)程一致，進(jìn)而導(dǎo)致整個(gè)水池區(qū)域的平均流量和湍流強(qiáng)度出現(xiàn)一些變化。

2）對(duì)于致動(dòng)盤(pán)水輪機(jī)模型，在TSR=5.5 時(shí)，y=0 m，z=1 m 截面上速度與湍流強(qiáng)度實(shí)驗(yàn)值與仿真值兩者誤差分別為3%到5%。隨著尾流的擴(kuò)大，x=0 m 截面時(shí)中心速度最大，x=4.0 m 截面時(shí)，其中心速度最小。且由于湍流和機(jī)架作用尾流中心處的速度逐漸下降后，兩側(cè)尾流速度逐漸恢復(fù)。

應(yīng)用推力T作為致動(dòng)盤(pán)模型附加源項(xiàng)，使仿真模擬具有計(jì)算時(shí)間少所需計(jì)算資源少優(yōu)點(diǎn)，并且也能精準(zhǔn)捕捉尾流場(chǎng)變化。圓形拖曳池相較于方形拖曳池在造波性能上有范圍廣、方向多、符合真實(shí)環(huán)境等優(yōu)點(diǎn)，但有占用面積大、建設(shè)費(fèi)用較高、不利于大型實(shí)驗(yàn)等缺點(diǎn)。由于圓形拖曳池中心區(qū)較小，不利于開(kāi)展大規(guī)模陣列水輪機(jī)實(shí)驗(yàn)，因此后續(xù)將利用方形拖曳池對(duì)水平軸水輪陣列間距機(jī)進(jìn)行實(shí)驗(yàn)研究。