劉葳興，劉 磊，陳雨龍，王 鑫，張之陽

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

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

1 實驗裝置與方法

1.1 試驗拖曳池

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

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

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

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

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

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

2 致動盤模型

致動盤方法基于一維動量理論，將水輪機葉輪視為一個可穿透的等效面積圓盤區(qū)域代替[14]。同時，將水輪機致動盤周圍流場視為穩(wěn)態(tài)、一維、不可壓流場，并以體積力的形式實現(xiàn)致動盤對流場的作用[15]。

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

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

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

Δx為致動盤厚度，V為致動盤體積。其動量方程中，致動盤體積力源項表示為[17]

本研究采用RNGk-ε湍流模型，k代表湍流動能，ε代表這個能量的耗散率。

3 數(shù)值模擬

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

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

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

圖2 渦輪機網(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ì)量流出或無流流動。將圖4 拖曳池中的28 個葉輪單元配置為入口、出口及壁面，質(zhì)量流量大小根據(jù)表2 設(shè)定。表2 中質(zhì)量流量（X軸流向）正值為流入，負(fù)值為流出。流入之和與流出之和為相同數(shù)值，則保持質(zhì)量守恒。水池壁面速度設(shè)置為0 m/s，壁面函數(shù)用于k，ε。流體域頂部設(shè)置為無滑移壁面。除速度外，其余初始條件皆設(shè)置為邊界條件，初始速度設(shè)置為0 m/s，動力粘度ν設(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é)果分為兩部分，一是拖曳池數(shù)值仿真，即在無水輪機情況下水流在拖曳池中自由發(fā)展；二是在致動盤模型下進行仿真。將兩者仿真結(jié)果對比，分析水輪機尾流虧損情況，并為后續(xù)研究做基礎(chǔ)。

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

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

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

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

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

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

4.2 致動盤模型仿真

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

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

其中，T為推力，ρ流體密度,U為來流速度。

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

其中，P為水輪機實際功率，A為致動盤面積。

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

圖8 水輪機仿真Fig.8 Tidal turbine simulation contour

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

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

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

5 結(jié)論

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

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

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

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