何 聰,史宏達(dá),張 嶔,趙昌宇,翟仁彬

(中國(guó)海洋大學(xué) 工程學(xué)院,山東 青島266100)

潮流能因可預(yù)測(cè)性高、儲(chǔ)量豐富而被廣泛研究。水平軸潮流能水輪機(jī)能很好地轉(zhuǎn)換利用潮流能,但高昂的發(fā)電成本影響了潮流能的商業(yè)發(fā)展。采用多機(jī)組潮流能水輪機(jī)陣列能有效節(jié)省發(fā)電成本、增加供電量,通過對(duì)尾流結(jié)構(gòu)研究可以確定相鄰水輪機(jī)之間的最佳距離,從而獲得更高的轉(zhuǎn)換效率[1]。同時(shí),水輪機(jī)尾流對(duì)海床沖刷、積淤等也有重要影響[2]。因此研究潮流能水輪機(jī)的尾流結(jié)構(gòu)有重要意義。

在潮流能尾流場(chǎng)研究中,Gaurier等[3]利用實(shí)驗(yàn)對(duì)比了不同湍流下的水輪機(jī)性能,結(jié)果表明低流速下水輪機(jī)荷載會(huì)受內(nèi)部湍流結(jié)構(gòu)的影響。張玉全等[4]討論了水輪機(jī)安裝高度對(duì)尾流擴(kuò)散的影響,提出增加水輪機(jī)與池底的距離可以加快尾流恢復(fù)速度,減小尾流長(zhǎng)度。但尾流復(fù)雜的物理規(guī)律沒有被完全揭示,隨著科學(xué)技術(shù)進(jìn)步,數(shù)值模擬因成本低、效率高逐漸應(yīng)用于潮流能水輪機(jī)尾流的空間演變過程研究。Leroux等[5]使用RANS-SST模型對(duì)潮流能水輪機(jī)進(jìn)行了數(shù)值模擬研究,結(jié)果表明穩(wěn)態(tài)和瞬態(tài)的方法在潮流能水輪機(jī)性能模擬方面精確度相似,都能較好預(yù)測(cè)推力系數(shù)和功率系數(shù),但是在對(duì)尾流的數(shù)值模擬中,瞬態(tài)模擬的結(jié)果更加符合實(shí)驗(yàn)結(jié)果。Afgan等[6]對(duì)比了基于RANS和LES模型對(duì)潮流能水輪機(jī)尾流場(chǎng)的數(shù)值模擬差異,結(jié)果表明LES模型可以捕捉到葉片尖端旋渦及解釋支撐塔架對(duì)尾流的相互影響,而RANS未捕捉到這2種旋渦。相比RANS模型,LES模型對(duì)流場(chǎng)的描述更加準(zhǔn)確詳細(xì)。但需要昂貴的計(jì)算成本,且受限于計(jì)算設(shè)備。為了減少求解邊界層湍流結(jié)構(gòu)的計(jì)算量,DES模型被提出并后續(xù)優(yōu)化為DDES模型。

本文基于DDES方法對(duì)潮流能水輪機(jī)尾流場(chǎng)進(jìn)行數(shù)值模擬,通過與實(shí)驗(yàn)結(jié)果對(duì)比,驗(yàn)證數(shù)值模擬的計(jì)算精度和可靠性,敘述數(shù)值模擬的網(wǎng)格類型、模型參數(shù)。重點(diǎn)分析水輪機(jī)尾流場(chǎng)的典型物理量的分布特征,以及研究水輪機(jī)的渦結(jié)構(gòu)發(fā)展過程。

1 幾何模型與計(jì)算區(qū)域

潮流能水輪機(jī)幾何模型采用Payne等[7]設(shè)計(jì)的三葉片水平軸水輪機(jī)模型(圖1)。轉(zhuǎn)盤半徑(R)為0.6 m,葉片翼型為NACA 63-8XX。水輪機(jī)尾艙由兩部分組成,直徑(D)由0.12 m增加至0.16 m。數(shù)值模擬的計(jì)算域大小與實(shí)驗(yàn)水槽尺寸相符,長(zhǎng)11D,寬3.3D,高1.7D(圖2)。計(jì)算域原點(diǎn)位于水輪機(jī)轉(zhuǎn)子頂部中心。

圖1 水輪機(jī)模型Fig.1 Model of the water turbine

圖2 計(jì)算域尺寸Fig.2 Di mension of the computation domain

對(duì)入口的邊界條件設(shè)定為速度入口,(x,y,z)向流速分別為(0.81,0,0)m·s-1,湍流強(qiáng)度為12%。頂部、兩側(cè)和底部的邊界類型為對(duì)稱邊界條件“sy mmetry”,出口處設(shè)置為壓力出口,水輪機(jī)的葉片和尾艙部分設(shè)置為無(wú)滑移壁面邊界類型。

2 數(shù)值模型建立

2.1 湍流模型

DDES模型較DES模型[8]引入一個(gè)延遲函數(shù),重新構(gòu)造了DDES的長(zhǎng)度尺度,同時(shí)考慮了網(wǎng)格尺度和渦黏場(chǎng),避免了從RANS到LES的切換太過靠近壁面,也防止模型預(yù)測(cè)過早分離。本文基于SST模型的DDES-SST方法[9],其控制方程為

式中:ρ為流體密度;U為速度;k為瞬時(shí)湍流動(dòng)能;lDDES為DDES長(zhǎng)度尺度;ω為比耗散率;σk=0.85、σω=0.5、σω2=0.856、a1=0.31、α=1、β=0.09,均為模型系數(shù);S為應(yīng)變力張量;Pk、Pw為結(jié)果項(xiàng);P為壓力;μt為湍流渦黏度;F1、F2為SST的混合函數(shù)。式(1)中DDES的長(zhǎng)度尺度lDDES為

式中:ILES、IRANS分別為L(zhǎng)ES、RANS的長(zhǎng)度尺度。

DDES的長(zhǎng)度尺度l 和渦黏性場(chǎng)、時(shí)間有關(guān),其作用能使RANS湍流模型具有自我延續(xù)功能。如果數(shù)學(xué)函數(shù)fd表明某點(diǎn)位于邊界層內(nèi),則使用長(zhǎng)度尺度lDDES的模型可以判斷轉(zhuǎn)變?yōu)長(zhǎng)ES模式的時(shí)機(jī)。當(dāng)出現(xiàn)大范圍的分離流動(dòng)時(shí),fd值就會(huì)從0轉(zhuǎn)變?yōu)榈絃ES模型。流體發(fā)生分離后從RANS到LES的轉(zhuǎn)換比DES方法更迅速。DES結(jié)合了RANS和LES的優(yōu)勢(shì)同時(shí)相比LES計(jì)算量減少,而DDES改善了DES存在的灰色區(qū)域和?；瘧?yīng)力不足的問題。

2.2 網(wǎng)格劃分

結(jié)構(gòu)化網(wǎng)格能夠?qū)ξ擦鲌?chǎng)信息捕捉效果更好,網(wǎng)格較為規(guī)整。非結(jié)構(gòu)化網(wǎng)格對(duì)復(fù)雜區(qū)域處理更為高效,同時(shí)也能保證數(shù)值模擬的精確性[10]。因此本文利用Point wise18網(wǎng)格劃分軟件,采用混合網(wǎng)格類型對(duì)潮流能水輪機(jī)進(jìn)行網(wǎng)格生成。即水輪機(jī)葉片面網(wǎng)格、葉片周圍的體網(wǎng)格以及水輪機(jī)尾流區(qū)域采用結(jié)構(gòu)化網(wǎng)格,轉(zhuǎn)子頂部、尾艙表面網(wǎng)格使用非結(jié)構(gòu)化三角形網(wǎng)格(圖3)。借助軟件中的“Voexl”算法完成計(jì)算域中體網(wǎng)格的填充,該算法以四面體、金字塔、棱柱等類型填充過渡區(qū)域,在尾流區(qū)域填充六面體結(jié)構(gòu)化網(wǎng)格(圖4)。網(wǎng)格拓?fù)浣Y(jié)構(gòu)分為動(dòng)域與靜域兩部分,動(dòng)域包含水輪機(jī)葉片與轉(zhuǎn)子,靜域?yàn)榱鲌?chǎng)其余部分(圖5)。此拓?fù)浣Y(jié)構(gòu)有利于水輪機(jī)進(jìn)行定常的多重參考系、非定常的滑移網(wǎng)格方法計(jì)算。

圖3 水輪機(jī)面網(wǎng)格劃分Fig.3 Surface grid divisions of the water turbine

圖4 Voexl算法對(duì)葉片附近填充的體網(wǎng)格Fig.4 The volu me grids filled near the blades based on the Voexl algorith m

圖5 計(jì)算域網(wǎng)格Fig.5 Grids of the co mputation do main

2.3 網(wǎng)格無(wú)關(guān)性分析

以粗網(wǎng)格、中網(wǎng)格、細(xì)網(wǎng)格類型,設(shè)計(jì)了3種密度的水輪機(jī)網(wǎng)格數(shù)量,詳見表1。在加密動(dòng)域的同時(shí),靜域網(wǎng)格數(shù)也在增加,使得動(dòng)域、靜域交界面網(wǎng)格尺寸大小相近,以保證網(wǎng)格連續(xù)性和良好的過渡性,提高數(shù)值模擬計(jì)算的精確性。

表1 5種網(wǎng)格的靜域、動(dòng)域網(wǎng)格數(shù)Table 1 The grid nu mber of the static and the dyna mic do mains of 5 kinds of grids

采用湍流“Spalart-All maras”模型進(jìn)行穩(wěn)態(tài)計(jì)算,停止標(biāo)準(zhǔn)為迭代5 000次,在迭代3 000次時(shí)水輪機(jī)推力值趨于穩(wěn)定。水輪機(jī)網(wǎng)格在計(jì)算得出壁面的最大Y+值在5以內(nèi),出現(xiàn)在葉尖處,最小首層網(wǎng)格厚度為1×10-7m。隨著水輪機(jī)網(wǎng)格的不斷加密,推力系數(shù)也在增加。當(dāng)網(wǎng)格加密至4 182 W時(shí),推力系數(shù)監(jiān)測(cè)值變化幅度減小,而粗網(wǎng)格和細(xì)網(wǎng)格的推力差距較大。綜合計(jì)算精度、計(jì)算資源成本,本文最終選擇細(xì)網(wǎng)格作為數(shù)值模擬的最終網(wǎng)格。

3 計(jì)算結(jié)果與討論

3.1 水輪機(jī)性能結(jié)果驗(yàn)證

水輪機(jī)性能與流場(chǎng)參數(shù)定義如下:

式中:RTS為轉(zhuǎn)速比;CT,CP分別為水輪機(jī)的推力系數(shù)和功率系數(shù);P為功率;ETK為湍流動(dòng)能(m·s-1)2;ω為水輪機(jī)旋轉(zhuǎn)速度(rad·s-1);R為水輪機(jī)半徑(m);U0為入口流速(m·s-1);T為轉(zhuǎn)子推力(N);ρ為水的密度(kg·m-3);A為轉(zhuǎn)盤面積(m2);Q為水輪機(jī)轉(zhuǎn)矩(N·s);Ux、Uy、Uz分別為x、y、z流向的瞬時(shí)速度與平均速度差值(m·s-1)。

以細(xì)網(wǎng)格進(jìn)行定常計(jì)算,設(shè)置水輪機(jī)RTS分別為5.37,5.52,5.75,5.98,6.52,6.97和7.51。通過功率系數(shù)、推力系數(shù)式(6)～式(7)計(jì)算得出水輪機(jī)在不同轉(zhuǎn)速下的推力系數(shù)與功率系數(shù)。由圖6可見,數(shù)值模擬結(jié)果低于實(shí)驗(yàn)結(jié)果的7%,但是仍處于實(shí)驗(yàn)采集數(shù)據(jù)的相對(duì)誤差10%范圍內(nèi)。由于本文主要對(duì)尾流場(chǎng)發(fā)展規(guī)律進(jìn)行研究,未考慮實(shí)驗(yàn)中的支撐結(jié)構(gòu),該部分誤差可能來(lái)源于支撐結(jié)構(gòu)的影響[11]。

圖6 實(shí)驗(yàn)與數(shù)值模擬的功率系數(shù)和推力系數(shù)對(duì)比Fig.6 Co mparison of power coefficient and thrust coefficient bet ween the experi ments and the nu merical si mulations

3.2 水輪機(jī)尾流結(jié)果驗(yàn)證

在入口x向流速為0.81 m·s-1、轉(zhuǎn)速為9.01 rad·s-1工況時(shí),非定常模型計(jì)算的時(shí)間步長(zhǎng)為9.69×10-4s,每步最大迭代次數(shù)8次。達(dá)到停止時(shí)間32 s后,水輪機(jī)在旋轉(zhuǎn)10圈后流場(chǎng)達(dá)到穩(wěn)定狀態(tài)。取數(shù)值模擬結(jié)果中水輪機(jī)旋轉(zhuǎn)最后5圈數(shù)據(jù)進(jìn)行相位平均處理。由此得出水輪機(jī)尾流場(chǎng)中X=0.1D,0.5D和1.0D截面上監(jiān)測(cè)點(diǎn)的平均速度,并與實(shí)驗(yàn)監(jiān)測(cè)數(shù)據(jù)[7]進(jìn)行對(duì)比。在Y=0處,因?yàn)槲才摰拇嬖?X=0.1D,0.5D截面下,尾艙直徑的部分實(shí)驗(yàn)與數(shù)值數(shù)據(jù)均無(wú)法采集。由圖7可見,數(shù)值模擬的尾流速度與實(shí)驗(yàn)尾流的截面數(shù)據(jù)吻合良好,下游數(shù)值結(jié)果中出現(xiàn)的差異主要與數(shù)值模擬中使用的簡(jiǎn)化有關(guān),其中沒有對(duì)支撐結(jié)構(gòu)進(jìn)行建模,而尾艙被短圓柱代替。上述結(jié)果說(shuō)明數(shù)值計(jì)算模型的計(jì)算精度較高,對(duì)后續(xù)尾流場(chǎng)及渦特性分析具備一定的可信度。

圖7 3個(gè)截面下數(shù)值模擬與實(shí)驗(yàn)尾流監(jiān)測(cè)點(diǎn)速度對(duì)比Fig.7 Comparison of the wake speeds at the monitoring points of three sections bet ween the experi ments and the numerical si mulations

3.3 水輪機(jī)尾流特性分析

本文以不同轉(zhuǎn)速、不同流速設(shè)計(jì)了4種工況,以此進(jìn)行尾流場(chǎng)結(jié)果的對(duì)比分析,并研究尾流場(chǎng)的發(fā)展規(guī)律。同時(shí)借助高性能計(jì)算集群,以192核的計(jì)算資源對(duì)潮流能水輪機(jī)單個(gè)工況并行計(jì)算,各算例詳細(xì)信息及計(jì)算所耗時(shí)間見表2。

表2 數(shù)值模擬中各工況詳細(xì)信息及計(jì)算耗時(shí)Table 2 The detailed infor mation and co mputation ti me of each working conditions in the nu merical si mulations

湍流是非常復(fù)雜的隨機(jī)過程,對(duì)湍流動(dòng)能的分析可以了解流場(chǎng)湍流隨時(shí)間和空間變化的情況。為貼近真實(shí)的海洋流場(chǎng)情況,入口邊界至X=3D處的湍流動(dòng)能值較高。對(duì)湍流動(dòng)能數(shù)值高于10-4的區(qū)域進(jìn)行遮擋,以便對(duì)尾流場(chǎng)湍流動(dòng)能觀察更為清晰。由圖8可見,4種工況中流場(chǎng)湍流動(dòng)能均沿Z=0對(duì)稱分布,低流速流場(chǎng)在X=5D處湍流動(dòng)能開始減弱。此后,流場(chǎng)的發(fā)展以自由剪切流與尾流混合發(fā)展為主。湍流動(dòng)能在X=8.5D處向中心聚集。高流速流場(chǎng)中,尾艙后的湍流動(dòng)能數(shù)值為低流速流場(chǎng)3倍,高轉(zhuǎn)速下的流場(chǎng),湍流動(dòng)能在發(fā)展至X=6D時(shí),湍流動(dòng)能有增加的趨勢(shì),在X=9.5D處湍流動(dòng)能明顯增加。這說(shuō)明在高轉(zhuǎn)速情況下,流場(chǎng)與水輪機(jī)轉(zhuǎn)速后的尾流與自由流存在混合現(xiàn)象,增強(qiáng)了尾流的湍流效應(yīng),從而導(dǎo)致湍流動(dòng)能增加。低轉(zhuǎn)速下,湍流動(dòng)能從X=5.5D開始減弱,但流場(chǎng)發(fā)展仍以水輪機(jī)旋轉(zhuǎn)產(chǎn)生的尾流為主。

圖8 4種工況下Y=0截面的湍流動(dòng)能云圖Fig.8 Cloud map of t he tur bulent kinetic energy at Y=0 section under the four wor king conditions

圖9 中給出了4種工況下的渦量云圖。渦量是描述流體的旋轉(zhuǎn)運(yùn)動(dòng)特性的物理量,能表征渦特性及渦分離、脫落演化過程,由對(duì)速度的旋度計(jì)算得來(lái)。x,y,z方向渦量的計(jì)算公式為

式中:u、v、w分別為x、y、z方向的速度。渦量(單位為S-1)值大小的計(jì)算方法為‖ω‖=

由圖9可見,從Y=0截面觀察到水輪機(jī)葉尖后產(chǎn)生2個(gè)主渦帶,轉(zhuǎn)軸后為2個(gè)次渦帶,輪轂渦在經(jīng)過尾艙后發(fā)展成為一個(gè)渦帶。由于葉片的高速旋轉(zhuǎn),而葉尖處形成葉尖渦,并逐步脫落形成葉尖脫落渦。4種工況下的葉尖脫落渦以規(guī)則圓形向后發(fā)展。并且沿流場(chǎng)方向直徑增加,渦結(jié)構(gòu)失穩(wěn)直至消散。在工況1與工況2中能清晰觀察到繞尾艙旋轉(zhuǎn)的輪轂渦。由水輪機(jī)葉片根部旋轉(zhuǎn)與尾艙的相互作用產(chǎn)生,并從尾艙末端脫落。高流速流場(chǎng)葉尖脫落渦在X=7D處開始減弱、消散,而低流速流場(chǎng)葉尖脫落渦在X=5D處開始消散。

圖9 4種工況下的渦量云圖Fig.9 Cloud map of the vorticity under the f our working conditions

圖10 為4種工況下的尾流場(chǎng)渦量等值面圖,通過對(duì)水輪機(jī)的三維渦量等值面圖分析,可以觀察到水輪機(jī)中渦結(jié)構(gòu)的詳細(xì)發(fā)展過程。為了更清晰地觀察流場(chǎng)渦結(jié)構(gòu)的變化,以不同渦量值形成等值面,設(shè)置工況1～工況4的顯示渦量值分別為3.0、3.8、3.5和5.0。圖中展現(xiàn)的葉尖脫落渦最終都是以橢圓形狀脫落,但是脫落過程不同,低流速流場(chǎng)較高流速流場(chǎng)的脫落過程往往更加迅速,葉尖脫落渦形成的渦環(huán)連續(xù)性更強(qiáng),渦環(huán)之間的距離減少0.1D。高流速流場(chǎng)中,在葉片旋轉(zhuǎn)與流場(chǎng)作用下,轉(zhuǎn)盤中間形成的渦數(shù)量更多。輪轂渦在低流速流場(chǎng)下,繞尾艙旋轉(zhuǎn)發(fā)展的跡象更突出,能觀察到輪轂渦旋轉(zhuǎn)時(shí)形成的渦環(huán)。輪轂渦在經(jīng)過直徑變大的尾艙時(shí),部分輪轂渦發(fā)生破碎。流場(chǎng)速度增加時(shí),葉片與尾艙、流場(chǎng)之間的相互作用增強(qiáng),輪轂渦的發(fā)展情況更為復(fù)雜,渦破碎現(xiàn)象更為明顯。當(dāng)流速一定時(shí),轉(zhuǎn)速增加,葉尖渦脫落形成的渦環(huán)直徑增加。具體變化為工況1渦環(huán)直徑由0.18D增加至與工況3的0.2D,工況2渦環(huán)直徑由0.19D增加至與工況4的0.22D。

4 結(jié) 論

本文采用DDES方法對(duì)水輪機(jī)尾流場(chǎng)進(jìn)行數(shù)值模擬研究,通過對(duì)比不同入流與水輪機(jī)轉(zhuǎn)速尾流場(chǎng)的數(shù)值模擬結(jié)果,得出以下結(jié)論:

①4種工況下湍流動(dòng)能均沿Z=0軸對(duì)稱分布,轉(zhuǎn)速一定時(shí),流速越大湍流動(dòng)能越高。流速一定時(shí),轉(zhuǎn)速越大流場(chǎng)的自由流與尾流混合現(xiàn)象越明顯。

②水輪機(jī)尾流場(chǎng)渦結(jié)構(gòu)主要由葉尖渦、葉尖脫落渦和輪轂渦組成。流場(chǎng)的入流速度與轉(zhuǎn)速會(huì)影響尾流場(chǎng)中渦結(jié)構(gòu)發(fā)展距離與渦量形成大小。轉(zhuǎn)速一定時(shí),流場(chǎng)流速越大,葉尖脫落渦、輪轂渦的發(fā)展距離越遠(yuǎn)。流速一定時(shí),轉(zhuǎn)速越大,渦的發(fā)展距離越短。

本文對(duì)比了不同工況潮流能水輪機(jī)尾流渦結(jié)構(gòu)發(fā)展變化規(guī)律,研究結(jié)論可為潮流能水輪機(jī)的陣列布局提供技術(shù)支持,對(duì)未來(lái)環(huán)境保護(hù)、海床沖刷研究奠定基礎(chǔ)。同時(shí),研究結(jié)論可為潮流能水輪機(jī)的陣列布局提供技術(shù)支持,對(duì)未來(lái)環(huán)境保護(hù)、海床沖刷研究奠定基礎(chǔ)。