袁 鵬， 劉小棟， 司先才， 王樹杰， 譚俊哲， 鄭志爽

(1.中國海洋大學工程學院，山東 青島 266100； 2.青島市海洋可再生能源重點實驗室，山東 青島 266100)

隨著全球變暖等一系列由于化石燃料燃燒導致的環(huán)境問題以來，人們對可再生能源的開發(fā)和利用已經(jīng)愈加重視。其中，潮流能作為一種開發(fā)潛力巨大的可再生能源，在近些年已經(jīng)逐步引起了人們的關注。英國[1]、美國[2]、法國[3]等國家都已開始了潮流能的研究，在2009年歐洲海洋能中心(EMEC)頒布了第一份潮流能開發(fā)標準體系[4]并建立了較大規(guī)模的潮流能設備海試試驗基地。到目前為止，各類潮流能水輪機也在不斷地研究和測試中，其發(fā)展水平已接近產(chǎn)業(yè)化[5]。

相比于國外而言，我國潮流能技術也同樣發(fā)展迅速[6]。近些年，在國內(nèi)的許多水域進行了較多潮流能資源評估的研究，并對潮流能轉化裝置進行了研制[7]。在開發(fā)潮流能資源時，潮流能轉化裝置的安裝可能會對其周圍的自然水流模式產(chǎn)生影響，從而改變沉積物的分布和輸運，進而影響物理環(huán)境，尤其是大型商用陣列的布放可能逐漸改變海洋物理環(huán)境，因此在進行大規(guī)模潮流能開發(fā)時有必要對潮流能轉化裝置陣列對水域潮流特性的影響進行研究。

以青島齋堂島及附近水域為研究對象，該水域的潮汐類型屬于正規(guī)半日潮，最大潮差達到4.6 m，平均差約為2.8 m。其中，在齋堂島東南部水域水深為35～40 m且水底相對較為平坦，大潮期間峰值流速能超過1.9 m/s，是我國北方比較有代表性的潮流能資源區(qū)。近些年來，許多關于潮流能的研究都是圍繞該水域進行的：在2012年，李華軍等[8]通過Flux方法對其瑯琊臺海峽的潮流能資源進行了評估；同年，邵萌等[9]對在該水域建立500 kW海洋能獨立電力系統(tǒng)示范工程進行了規(guī)劃；在2013年，史廉博等[10]基于FVCOM建立三維水力學模型并對水域潮流特性進行分析；在2014年，紀合盼等[11]通過實測水文數(shù)據(jù)對潮流能發(fā)電站的建設進行初步規(guī)劃；在2015年，陳婭玲等[12]采用動量損失的方法模擬水輪機，并對該水域的水輪機陣列方式進行了研究；在2017年，林杰等[13]基于Delft 3D研究了該水域的漩渦等因素對潮流能提取可能帶來的影響。從現(xiàn)有文獻來看，目前的研究多集中于對裝置發(fā)電能力的優(yōu)化以及區(qū)域潮流能資源的評估，而對潮流能裝置對海洋環(huán)境的影響研究還很少。

本文基于Delft 3D-Flow建立齋堂島及附近水域的三維水力學模型，并通過實測數(shù)據(jù)對模型進行驗證。經(jīng)過驗證后的模型，在合理確定該水域中水輪機陣列放置位置[14]的基礎上，利用等效[15]方式模擬水平軸潮流能水輪機陣列，對加入水輪機陣列后的水力學模型進行計算，分析水輪機陣列對區(qū)域流場的水位和流速等影響，為未來在此水域開發(fā)潮流能資源提供參考。

1 潮流場水力學模型

1.1 建立模型

本次研究采用Delft 3D的Flow模塊建立水力學模型。基于淺水方程和Boussinesq假設，該模型采用有限差分法求解斜壓下的N-S輸運方程，既可以用于二維模型，也可以用于三維模型計算[16]。其主要基于三組控制方程，分別是：

質量守恒方程：

(1)

動量守恒方程：

(2)

(3)

輸運方程：

(4)

其中:ζ表示水位(總水深)；d表示到參考平面的當前水深(凈水深)；U、V分別代表在x，y方向上的速度分量；Q代表單位面積上質量源強度；f為科氏力參數(shù)；vh是動態(tài)水平渦流粘度；ρ0是參考密度；ρ′是不規(guī)則密度；τsx、τsy是作用在海平面上的風壓分量；τbx、τby為底部的剪應力分量。在輸運方程中，c代表鹽度或者溫度；Dh是水平渦流消散度；λd代表一階衰減過程；g為重力加速度；R為單位面積上的源項。

本次研究的模擬水域為青島齋堂島及其附近水域。齋堂島隸屬于山東省青島市，位于瑯琊臺東南方向海中，地理坐標為北緯35°38′，東經(jīng)119°55′(見圖1)。

為提高模型輸出結果的計算精度和減少計算機資源利用，在建模過程采用網(wǎng)格嵌套(Nesting)方法，分別建立2個水域范圍和網(wǎng)格分辨率不同的模型(見圖2，圖中黑色區(qū)域為大范圍水域網(wǎng)格范圍，紅色區(qū)域為小范圍網(wǎng)格范圍，藍線為陸線)。

圖1 齋堂島地理位置

大范圍水域所建立的模型由其開邊界上通過OTIS提供的8個潮汐調(diào)和常數(shù)(M2、2、N2、K1、K2、O1、P1、Q1)進行驅動，考慮到OTIS的輸出的潮汐調(diào)和常數(shù)的分辨率(1/30(°))，需要建立足夠大小的水域計算網(wǎng)格范圍保證足夠數(shù)量的調(diào)和常數(shù)作為邊界條件，從而確保模型仿真結果的可靠性。網(wǎng)格分辨率由水平方向上的水深數(shù)據(jù)分辨率決定，網(wǎng)格尺寸選取600 m×600 m并在齋堂島區(qū)域局部加密到120 m×120 m。

小范圍水域模型的建模水域范圍由本次研究所需的研究水域范圍決定，其驅動使用大范圍水域模型輸出的水位時間序列作為輸入。網(wǎng)格尺寸選取36 m×36 m并局部加密到18 m×36 m以保證后續(xù)對水輪機陣列等效模擬的精度。

圖2 模型網(wǎng)格

計算時大范圍水域采用二維水力學模型，小范圍水域采用三維水力學模型。其中，小范圍的水力學模型在垂直方向上將水流分為10 sigma層(每層厚度為10%的水深)并采用k-ε湍流模型進行計算。在確保滿足Courant條件的前提下，大小模型的時間步分別設置為1和0.1 min。利用所建立的水力學模型模擬一個半月潮的潮汐周期并且不考慮風、浪等影響，模擬時間選擇為17 d，其中前2 d作為模型由靜水狀態(tài)開始運行的啟動時間。

1.2 潮流場水力學模型驗證

為確保計算結果的有效性，分別通過實測的水位數(shù)據(jù)對大范圍水域模型進行驗證和實測流速流向數(shù)據(jù)對小范圍水域模型進行驗證。

1.2.1 大范圍模型驗證 在大范圍水域模型中，選取日照港、董家口港和齋堂島潮位站3個站點2013—2014年的部分實測水位數(shù)據(jù)進行驗證。其中大范圍模型的驗證點坐標以及部分驗證結果分別見表1和圖3。

通過實測水位數(shù)據(jù)與模型模擬結果進行對比：日照港水位驗證的平均水位誤差為5.58%，局部最大誤差為13.06%；齋堂島的平均水位誤差為8.00%，局部最大誤差為15.10%；董家口港平均水位誤差為9.94%，局部最大誤差為13.86%。

表1 潮位驗證點經(jīng)緯度

1.2.2 小范圍模型驗證 在大范圍水域模型經(jīng)驗證后，將其輸出的水位時間序列作為小范圍水域模型的邊界條件進行輸入， 然后將小范圍水域模型的輸出結果與實測流速流向數(shù)據(jù)進行對比，并在此基礎上調(diào)節(jié)小范圍水域模型的參數(shù)從而提高模型輸出結果的精確度。該水域:4個流速驗證點的坐標和位置分別見表2、圖4。

圖3 齋堂島潮位站驗證結果

流速驗證Velocity verification站點1Site 1站點3Site 3站點4Site 4站點5Site 5經(jīng)度Longitude/°E119°55′39.9″119°55′46.5″119°55′39.2″119°56′8.79″緯度Latitude/°N35°37′27.8″35°37′27.3″35°37′18.5″35°36′12.88″

根據(jù)Rahman A等的研究[17]，水力學模型的輸出流速受水底摩擦系數(shù)影響較大。而本次建立水力學模型的水域，因為缺乏齋堂島附近水域的精確水底摩擦數(shù)據(jù)，此次仿真采取統(tǒng)一的水底摩擦系數(shù)。在此，分別選用Manning=0.017、0.019和0.021對該區(qū)域流速驗證進行對比。圖5為小范圍水域模型的流速誤差驗證。

將模型的仿真結果和實測流速進行對比可知(見圖5)：1號站點流速驗證離散型較大且整體仿真流速偏大，可能由于其離島較近受到島嶼尾流所引起的漩渦影響較大[13]且水深變化劇烈，所以仿真精度較低；4號站點和5號站點處水深變化較為平緩且距離陸地相對較遠，流速驗證準確性較高；3號站點由于仍處于齋堂島尾流的影響范圍之內(nèi)，相對于4號和5號站點流速驗證準確性略低。

此外，根據(jù)圖5對比水底摩擦Manning系數(shù)不同時模型輸出各個站點的均方根誤差(RMSE)和相關系數(shù)(r)的數(shù)值，由此可以看出可知當Manning取0.019時，各站點的RMSE相對較小且r較大，因此流速驗證精度最高。因此，在后續(xù)的仿真模擬中，模型的Manning系數(shù)都選用0.019。

圖4 流速驗證點位置

圖5 各站點流速驗證

2 陣列對區(qū)域潮流特性改變研究

2.1 水輪機陣列放置位置選擇

潮流能水輪機陣列的位置通過TSE方法[14]來確定。TSE方法主要用于水深有限的水域中對水輪機位置進行選址，該方法很好地兼顧了區(qū)域潮流能能量強度和水深因素，并引入了潮流能可開發(fā)度(Tidal Stream Exploitability, TSE)因子來表征某區(qū)域中適合放置水輪機的程度，TSE因子的求解公式為：

(5)

式中：V0和h0分別為該水域的典型流速和水深，在這里分別取1 m/s和20 m；Vf和Ve分別為該水域漲落潮期間的流速；h為漲落潮期間的局部水深；ξ為補償函數(shù)，其計算公式為：

(6)

(7)

(8)

其中：h為漲落潮期間的水深；Δh為最大潮差；h1和h2分別為水輪機安裝水深補償?shù)淖畲笊舷孪?，在這里分別取2、5 m。

根據(jù)本次研究所建立的小范圍水域模型可得到該水域漲落潮階段的流速分布(見圖6)。仿真結果顯示，漲落潮階段在齋堂島東南部區(qū)域和其與大陸間的瑯琊臺海峽處流速較大。在大潮期間，漲潮流速能超過2.0 m/s，落潮流速能超過1.8 m/s。

圖6 大潮漲落潮階段流速分布

采用TSE方法可計算出該水域的TSE因子分布(見圖7)。在齋堂島東南部區(qū)域，TSE因子較大，比較適合放置規(guī)模較大的水輪機陣列；而在大陸與齋堂島間的海峽處，雖然漲落潮期間流速較大，但水深相對較小，不適合放置較大規(guī)模的水輪機陣列。因此，考慮到實際水深等因素，將水輪機陣列初步選定布放于齋堂島東南部水域，水輪機陣列位置如圖7黑線內(nèi)所示，其中陣列放置區(qū)域的面積為540 m×360 m。

(黑線內(nèi)部為水輪機陣列的放置區(qū)域。Black line is the placement area of the turbine array.)

圖7 齋堂島附近水域TSE因子分布圖
Fig.7 TSE index around waters of Zhaitang Island

2.2 水輪機陣列設置方式

在本次研究中，水輪機陣列中的水輪機都假設為水平軸潮流能水輪機。在流場模型中加入潮流能水輪機依研究的目的不同有不同的處理方法。其中，一些研究通過修改水力學模型中計算網(wǎng)格局部網(wǎng)格節(jié)點上的水底摩擦系數(shù)來等效水輪機陣列所帶來的水動力影響[18]，不過此方法多用于二維水力學模型[19]；另外一些研究通過在水力學模型中計算網(wǎng)格的局部節(jié)點上添加動量損失方程來仿真潮流能水輪機對水域產(chǎn)生的影響[12]；此外，也有一部分研究通過修改水力學模型中的局部湍流模型系數(shù)來模擬水輪機對潮流場帶來的改變[20]。

本次研究通過上述第二種方法，采用動量損失的方法。該方法的具體過程如下：

在水流通過潮流能水輪機時，會產(chǎn)生一定程度的動量損失，根據(jù)公式(2)和(3)，其動量方程分別可變?yōu)椋?/p>

(9)

(10)

其中，Mx和My分別為水力學模型在U和V兩個方向上的動量源項。

在Delft 3D中，可以通過UDF文件多孔盤(porousplate)來實現(xiàn)對其水力學模型中局部網(wǎng)格節(jié)點上動量方程的源項修改，以實現(xiàn)對潮流場中水輪機的模擬。該方法已通過實驗與Delft 3D中的仿真結果進行了一定驗證[15]并在近年的不少研究中得到了應用[19]。

圖9 多孔盤的在計算網(wǎng)格中布置方向

使用多孔盤仿真潮流能水輪機時，除了需要計算出水流通過潮流能水輪機時產(chǎn)生能量損耗的大小，還需要考慮在實際情況下的水輪機迎流方向與建立的水力學模型中計算網(wǎng)格方向不匹配的問題，因此需要在水輪機放置位置上的計算網(wǎng)格中的U和V兩個方向上同時設置多孔盤(見圖9)。多孔盤中參數(shù)設置的計算公式為：

(11)

其中，動量耗散參數(shù)Closs設置的計算公式，在U方向上為：

(12)

(13)

Au=∑At|cosθ| 。

(14)

同樣，在V方向上：

(15)

(16)

Av=∑At|sinθ| 。

(17)

其中：γu和γv分別代表在U和V方向上單位網(wǎng)格上的動量損耗比例；CT為水輪機的軸向力系數(shù)，本次仿真假設CT為定值并且選取CT=0.8[21]；Au和Av分別代表在U和V方向上單位網(wǎng)格上的水輪機葉片掃掠面積；At代表在單位網(wǎng)格上的水輪機的葉片上的掃掠面積；Δx和Δy分別代表單位網(wǎng)格節(jié)點上在U和V方向上的距離；Δz代表每sigma層的高度；n代表多孔盤所在的層數(shù)。

2.3 水輪機陣列對區(qū)域水位流速影響

本次研究，所選擇布放水輪機陣列的區(qū)域平均水深在35 m左右，根據(jù)EMEC的標準[4]，水平軸水輪機葉片最底端需要高出水底5 m以上，水輪機葉片最頂端需要低于該水域天文潮水位最低點5 m以上，因此選用葉片直徑為18 m的水輪機，其中水輪機轉子中心位置在距離海床高度15 m處。

在所選區(qū)域放置4排水輪機，每排間距離為10D(D為水輪機轉子直徑)；每列中相鄰水輪機轉子中心間距離分別選用1D、2D、3D、4D和5D，采用交錯布置的方式。在水輪機陣列布放面積確定的條件下，以上各陣列布置方式中可布放水輪機臺數(shù)分別為80、40、28、20和16臺。分別將未放置水輪機和放置不同排布方式的水輪機陣列的潮流場模擬結果進行對比可得到的半月潮大潮期間水位和流速矢量對比如圖10所示。

仿真結果顯示：漲落潮期間，水位會在水輪機陣列迎流方向上有一定程度的上升，且每經(jīng)過一排水輪機后會有一定的下降，該現(xiàn)象隨陣列中水輪機密度的增加而愈加顯著。其中，當一定面積陣列中水輪機密度達到最大時，其對區(qū)域水域的改變最大不超過3 cm，且隨其與陣列距離的增加影響幅度逐漸減小。另外，陣列對區(qū)域水位影響范圍較大，將會影響到島嶼與大陸間海峽處的水位(見圖10(a))。

相比對水位的影響，水輪機陣列的布放密度對流速影響十分明顯，其所帶來的最大流速下降在尾流影響區(qū)域能夠超過0.45 m/s，且影響距離超過2 km(見圖10(a))。而當選擇較為稀疏的水輪機陣列時(行間距大于4D時)，其對潮流場的影響較小，流速改變不到0.1 m/s(見圖10(e))。此外，在該水域加入水輪機陣列后將會對由島嶼尾流所引起的漩渦[13]帶來一定改變，此種改變所帶來的影響需要在未來進一步研究和評估。

((a) 列間距為1D,水輪機臺數(shù)共80臺；(b) 列間距為2D,水輪機臺數(shù)共40臺；(c) 列間距為3D,水輪機臺數(shù)共28臺；(d) 列間距為4D,水輪機臺數(shù)共20臺(e) 列間距為5D,水輪機臺數(shù)共16臺。(a) 80 turbines with 1D horizontal spacing； (b) 40 turbines with 2D horizontal spacing； (c) 28 turbines with 3D horizontal spacing； (d) 20 turbines with 4D horizontal spacing； (e) 16 turbines with 5D horizontal spacing.)

圖10 不同陣列布置方式下的水輪機陣列對區(qū)域水位和流速影響

Fig.10 Effect of turbine array on water level and velocity in different layout

3 結語

本文以Delft 3D的Flow模塊為基礎，建立齋堂島及其附近水域的三維水力學模型。通過TSE方法確定了水輪機陣列的布放位置，并將水輪機對流場的影響等效成為動量損失的方式并組成不同排列方式的陣列加入到水力學模型中，從而對流場進行仿真并研究陣列對潮流場所產(chǎn)生的影響。

仿真結果顯示，在齋堂島東南部水域放置一定規(guī)模的水平軸潮流能水輪機陣列將會對該區(qū)域的流場帶來一定影響：陣列會對較大范圍水域的水位帶來一定改變，且影響程度隨陣列中水輪機密度的減少而降低，其中最大水位改變不超過3 cm；對區(qū)域流速的影響隨陣列中水輪機密度的增加而愈加顯著，當密度較大時陣列所引起的流速下降在其尾流區(qū)域能夠達到0.45 m/s，影響距離可延伸至其下游2 km處。