宋佳奇，王俊榮，2，王德志，康信龍

(1.中國海洋大學 工程學院，山東 青島 266100;2.山東省海洋工程重點實驗室，山東 青島 266100)

海洋風能資源豐富、風場穩(wěn)定，海上風能資源開發(fā)已成為未來新能源開發(fā)的熱點。Spar風電平臺具有水線面小、重心低等特點，因此擁有優(yōu)良的運動性能和穩(wěn)性是開發(fā)深海風能的重要裝備。

由于Spar風電平臺吃水較深，在海流作用下會導致其后方產(chǎn)生交替脫落的旋渦，從而形成順流向和橫流向的脈動壓力，使平臺產(chǎn)生大幅的水平運動，稱之為渦激運動(Vortex-Induced Motion，VIM)。Spar風電平臺一般為無側板單立柱平臺，該種平臺形式更易發(fā)生大幅渦激運動。Hywind Scotland Spar型6MW風電平臺是立柱直徑14.4 m、吃水78 m 的無側板單立柱Spar風電平臺，該類平臺在實際運行中易監(jiān)測到顯著的大幅渦激運動[1]。渦激運動對系泊系統(tǒng)的疲勞帶來挑戰(zhàn)，需對其運動現(xiàn)象及規(guī)律進行深入研究。

目前，海洋平臺渦激運動的主要研究方法有實尺度觀測、模型試驗和計算流體動力學(Computational Fluid Dynamics，CFD)數(shù)值計算等，其中計算流體動力學數(shù)值計算方法具有可操作性強、易重復和耗費少等優(yōu)點，已逐步成為渦激運動的主要研究方法。2005年Thiagarajan等[2]采用計算流體動力學數(shù)值計算方法計算對比了有、無側板圓柱的渦激運動頻域特點;2014年蘇云龍[3]對單根圓柱固定繞流進行了計算流體動力學數(shù)值計算，通過流場分析研究了旋渦脫落機理，但未考慮流固耦合;2015年趙偉文和萬德成[4]用計算流體動力學數(shù)值計算方法計算了Ur為5～9時有側板Spar平臺的渦激運動響應，并與模型試驗進行了比較;2016年Duan 等[5]采用模型試驗的方法研究了OC3 Spar型浮式風機的渦激運動特性，發(fā)現(xiàn)并分析了平臺橫蕩、縱蕩鎖頻現(xiàn)象;2018年翟佳偉等[6]研究了風浪流聯(lián)合作用下Spar基礎的運動特性，其渦激力由二維圓柱繞流計算流體動力學數(shù)值計算模型計算得到;2019年何佳偉等[7]計算模型尺度下的有、無側板Spar平臺的渦激運動過程，重點關注了螺旋側板的減渦效果。針對無螺旋側板Spar風電平臺，獲取平臺渦激運動響應幅值與來流約化速度關系曲線，是分析平臺渦激運動特性和評估系泊系統(tǒng)疲勞的基礎。

本文建立了無側板單立柱Spar風電平臺二維計算流體動力學數(shù)值計算模型，開展平臺在不同流速下的渦激運動數(shù)值計算，探究Re不相似對平臺渦激運動的影響，為平臺相關設計分析與試驗提供參考。

1 數(shù)值方法

1.1 控制方程

渦激運動是非定常、黏性流問題，本文假設流體為不可壓縮牛頓流體，采用雷諾平均法(Reynolds Average Navier-Stockes，RANS)處理納維-斯托克斯方程(N-S方程)，得到的控制方程包括:

1)連續(xù)方程

2)動量方程

1.2 湍流模型

本文所使用SSTk-ω湍流模型是一種雷諾平均法，最早由Menter等[8-9]提出。該模型通過剪切應力輸運(Shear Stress Transport，SST)過程將k-ω模型和k-ε模型兩湍流模型結合，在近壁面區(qū)域采用邊界條件處理能力良好的k-ω模型，在遠場采用自由剪切流處理較好的k-ε模型。SSTk-ω模型充分結合k-ω模型kε模型的優(yōu)點，提高了計算精度。模型中湍動能(k)和特定湍流耗散率(ω)方程分別為[8]:

式中:μt為湍流渦動力黏性;νt為湍流渦運動黏性;P為湍流生成項;為進行限制后的生成項;β*為修正系數(shù);β，γ，σk和σω均為封閉參數(shù);F1為混合函數(shù)，用以在近壁處激活標準k-ω模型，在遠場處激活k-ε模型;CDkw(cross-diffusion)為交叉擴散項。式(3)和式(4)右側前3項分別為生成項、耗散項和擴散項，式(4)右側第4項為交叉擴散項。其中:

式中:S ij為平均應力張量;δij為克羅內克函數(shù);CD*kω為修正后的交叉擴散項;Ω為旋度幅度;y為距最近壁面距離;a1為修正系數(shù);F2為混合函數(shù);封閉參數(shù)β、γ、σk和σω在近壁處和遠場處具有不同值，下標1為近壁處參數(shù)值，下標2為遠場處參數(shù)值，其中各參數(shù)值為0.5，β1=0.075、σk2=1.0、σω2=0.856、β2=0.028、β*=0.09、κ=0.41、a1=0.31。β、β*、γ、σk和σω在距離面不同位置具有不同的值，須加權平均得到式(3)和式(4)中的值，對任意參數(shù)α，有:

1.3 流固耦合

計算物體流固耦合運動需要采用動網(wǎng)格技術，對于平臺大幅度往復運動的問題，重疊網(wǎng)格方法可實現(xiàn)多個相互獨立網(wǎng)格之間產(chǎn)生無拘束相對運動，網(wǎng)格不變形，可穩(wěn)定地進行流固耦合計算。通過自編程序計算物體運動的動力積分，實現(xiàn)對重疊網(wǎng)格的運動控制和流固耦合。具體做法是:在每一時間步內，將計算流體動力學數(shù)值計算求解得到的流體激勵代入平臺運動微分方程求解物體運動響應，然后對重疊網(wǎng)格位置及流場進行更新，完成該時間步的流固耦合，流程見圖1。

圖1 渦激運動流固耦合求解流程Fig.1 Flow chart of solving the VIM fluid-solid coupling

1.4 渦激運動相關特征參數(shù)

本文重點關注如下浮式平臺渦激運動無因次參數(shù):

式(15)～式(20)中:St為斯特勞哈爾數(shù);A/D為無量綱振幅，其中D為結構特征長度(立柱直徑)，A為平臺運動幅值;CL為升力系數(shù)，CD為阻力系數(shù)，采用阻力系數(shù)均值和升力系數(shù)標準差(CL′)來表征平臺水動力系數(shù);ρ為流體密度;U為來流速度;fv為平臺渦脫頻率;fn為平臺靜水橫蕩頻率;Xmax和Xmin分別為平臺順流方向偏移最大值和最小值;F Y和F X分別為水動力在橫流向和順流向的分量;SL和SD分別為浸沒部分向升力和阻力方向的投影長度。

2 計算模型及工況

2.1 平臺模型

本文重點研究Spar風電平臺在海流作用下的渦激運動響應特征，不考慮空氣引起的渦激運動。實際尺度下，Spar風電平臺立柱直徑為10 m，吃水120 m，質量8 450 t。因渦激運動響應主要表現(xiàn)在水平面內運動，立柱截面不變，兼顧計算效率問題，故采用平臺二維渦激運動模型進行渦激運動計算。在計算流體動力學數(shù)值計算時采用1∶50的縮尺比，通過在計算流體動力學數(shù)值計算中定義流體黏度，可滿足雷諾數(shù)Re相似準則。

2.2 系泊系統(tǒng)

平臺模型采用等效線性水平系泊提供回復力，系泊系統(tǒng)由4根線性彈簧組成，布置方式如圖2所示，平臺立柱中心為坐標系原點，每根彈簧與坐標軸的夾角為45°，平臺橫蕩、縱蕩固有周期相等，實際尺度下，系泊剛度為2×104N/m，預張力為1.5×106N;系泊系統(tǒng)整體的阻尼比為0.06。

圖2 平臺系泊方式示意圖Fig.2 Sketch map of the mooring mode of the platform

2.3 計算區(qū)域、邊界條件與網(wǎng)格劃分

計算區(qū)域分為背景網(wǎng)格和前景網(wǎng)格兩部分(圖3)，背景網(wǎng)格為矩形區(qū)域，以平臺立柱直徑D為基準，矩形長36D，寬16D;前景網(wǎng)格為正方形區(qū)域，邊長7D。為保證有足夠空間使流場中旋渦充分發(fā)展，背景網(wǎng)格中平臺設置在靠近流體入口的一側，平臺中心距左側流體入口處8D。

圖3 計算區(qū)域幾何簡圖Fig.3 Geometry diagram of the computational domain

計算區(qū)域邊界條件設置左側為速度入口邊界;右側為壓力出口邊界，給定出口靜壓力p=0，對流動中的其他物理量均由流場內部值外差得到;上下邊界設置為自由滑移壁面;平臺立柱表面為無滑移壁面邊界，邊界上流體速度與壁面速度相同;前景網(wǎng)格邊界設置為重疊網(wǎng)格邊界。

對于湍流模型參數(shù)，k和ω的初始流場和入口邊界條件可以通過經(jīng)驗公式計算得到[9]，即

式中:I為湍流強度，νt/ν為湍流黏度與層流黏度之比。在無滑移壁面邊界處，k和ω的法向梯度為0，，其中y1為壁面第一層網(wǎng)格單元中心到壁面距離。

為研究網(wǎng)格收斂性，選取3種不同密度網(wǎng)格進行網(wǎng)格測試，網(wǎng)格設置見表1所示。當Re=1.42×104時，計算3套網(wǎng)格的和St，在保證計算結果穩(wěn)定的前提下，時間步長取Δt=0.02 s，結果如表1所示。網(wǎng)格2與網(wǎng)格3計算結果接近，繼續(xù)加密網(wǎng)格對精度提高的效果有限，故采用網(wǎng)格2進行后續(xù)數(shù)值計算，網(wǎng)格劃分如圖4所示。

圖4 計算區(qū)域整體網(wǎng)格劃分Fig.4 Grid division of the entire computational domain

表1 不同網(wǎng)格計算結果Table 1 Results calculated based on the different sets of grids

3 計算結果與分析

3.1 固定繞流計算

為了驗證數(shù)值計算方法的可靠性，首先對平臺固定繞流進行數(shù)值計算研究，本文選取4個工況，對平臺受力計算結果進行統(tǒng)計分析，工況設置及統(tǒng)計各工況結果所得的St及如表2。將本文計算所得的St同Blevins[10]總結的St-Re關系曲線、蘇云龍[3]數(shù)值計算的St結果對比(圖5)，本文及蘇云龍計算得到的St均略高于Blevins所得St-Re關系曲線;將本文計算所得的同Schlichting[11]所得關系曲線、蘇云龍[3]數(shù)值計算的結果對比(圖6)，其數(shù)值接近Schlichting所得關系曲線。本文數(shù)值模型計算結果處于合理范圍。

圖5 斯特勞哈爾數(shù)與雷諾數(shù)的關系Fig.5 Relationship between the Strouhal number and the Reynolds number

圖6 阻力系數(shù)均值與雷諾數(shù)的關系Fig.6 Relationship between the mean of drag coefficient and the Reynolds number

表2 平臺繞流工況及統(tǒng)計結果Table 2 Conditions and statistics of the flow around the Spar

3.2 渦激運動計算

本文探究平臺在Ur=2～14范圍內的渦激運動規(guī)律，計算模型與原型滿足Re相似準則，流體動力黏度μ=2.837×10-6Pa·s，此時平臺橫蕩和縱蕩的固有頻率fn=0.086 Hz，稱該平臺系統(tǒng)為系統(tǒng)1，具體工況設置如表3所示。

表3 系統(tǒng)1渦激運動計算工況Table 3 Conditions for VIM calculation of System 1

圖7為各工況平臺縱蕩、橫蕩運動時歷曲線及平臺運動軌跡圖。對比圖7中不同Ur下平臺的運動軌跡，得到平臺運動的變化特點和規(guī)律:Ur=2時，平臺運動軌跡呈“8”字形，但幅值較小;Ur=3 時，運動軌跡開始發(fā)散;Ur=4時，平臺運動軌跡總體呈現(xiàn)出“C”字形的運動軌跡，橫蕩、縱蕩幅值均較大;Ur=7時，平臺運動軌跡變?yōu)楠M長的“8”字形，橫蕩幅值較大，縱蕩幅值較小;Ur=10 時，平臺運動軌跡開始發(fā)散，運動軌跡變?yōu)樵谝欢ǚ秶鷥鹊耐鶑瓦\動;Ur＞10時，平臺運動范圍大幅縮小，橫蕩、縱蕩幅值均明顯減小。

圖7 不同約化速度下縱蕩、橫蕩運動時歷曲線及運動軌跡Fig.7 Time history and trajectory of the surge and sway motions at different reduced velocities

通過對平臺CL時歷曲線做傅里葉變換，得到各工況下旋渦脫落頻率(fv)。圖8為Spar風電平臺不同Ur下的橫蕩頻率(fsway)、渦激運動渦脫頻率(fv-VIM)和根據(jù)Blevins[10]總結的St計算得到的固定繞流渦脫頻率(fv-Blevins)對比圖。由圖8可見，在Ur=3.5～10范圍內，fv-VIM與fsway幾乎一致，且與fv-Blevins相差較大。fv-VIM出現(xiàn)了2次頻率鎖定現(xiàn)象:在Ur=3.5～4.5時，平臺發(fā)生了縱蕩共振，此時為縱蕩鎖頻，fv-VIM為0.048 Hz，縱向激勵力振蕩頻率與平臺縱蕩固有頻率相近;在Ur=5～10時，平臺發(fā)生了橫蕩共振，此時為橫蕩鎖頻，fv-VIM為0.085 Hz，橫向激勵力振蕩頻率與平臺橫蕩固有頻率相近。U r=10時fv-VIM和fsway開始出現(xiàn)差別，Ur＞10時fv-VIM明顯高于fsway，接近于fv-Blevins。

圖8 橫蕩頻率、渦激運動渦脫頻率和固定繞流渦脫頻率Fig.8 A contrast map of f sway，f v-VIM and f v-Blevins

對各工況結果進行統(tǒng)計分析，得到圖9所示的不同Ur下平臺縱蕩和橫蕩運動特征參數(shù)對比。圖9中X和Y為平臺縱蕩、橫蕩偏移平均值，平臺縱蕩、橫蕩運動幅值A X和A Y可由最大值減去均值近似計算得到。結合圖8和圖9可知:在縱蕩鎖頻時A X和Xmax顯著變大，Ur=4.5時A X達到最大;在橫蕩鎖頻時A Y/D＞1。鎖頻現(xiàn)象顯著影響平臺渦激運動幅值，平臺設計時應考慮縱蕩鎖頻對平臺運動響應的影響。

圖9 無量綱化平臺縱蕩、橫蕩運動的均值和最大值Fig.9 Nondimensionalized mean and maximum of the surge and sway motions at different reduced velocities

3.3 不滿足雷諾相似準則的影響

物理水池試驗無法滿足Re相似準則，實驗結果與原型運動結果會產(chǎn)生差異，本節(jié)通過計算流體動力學數(shù)值計算方法探究不滿足雷諾相似準則對渦激運動計算結果的影響。設置計算模型的流體動力黏度為1.003×10-3Pa·s，平臺系統(tǒng)橫蕩、縱蕩固有頻率為0.084 Hz，稱此平臺為系統(tǒng)2(表4)。

表4 系統(tǒng)2工況設置表Table 4 Condition settings for VIM calculation of System 2

圖10 2個系統(tǒng)阻力系數(shù)均值和升力系數(shù)標準差對比Fig.10 Comparison of and C′L between the two systems

圖11 相同約化速度下系統(tǒng)1和系統(tǒng)2平臺運動軌跡對比Fig.11 Comparison of Spar motion trajectory between System 1and System 2 at the same reduced velocity

4 結 論

本文采用計算流體動力學數(shù)值計算方法對一無側板單立柱Spar風電平臺渦激運動現(xiàn)象進行數(shù)值計算，研究平臺在不同Ur的渦激運動特征。計算采用SSTk-ω湍流模型，通過自編動力積分程序和重疊網(wǎng)格實現(xiàn)渦激運動的實時流固耦合，對平臺在Ur=2～14的渦激運動響應進行了計算，分析了渦激運動模式隨Ur變化的規(guī)律，定性分析了不滿足雷諾相似準則對渦激運動計算的影響。得出以下結論:

1)平臺在不同Ur下會經(jīng)歷2次鎖頻，在Ur=3.5～4.5時出現(xiàn)運動軌跡呈“C”字形的縱蕩鎖頻，在Ur=5～10出現(xiàn)運動軌跡呈“8”字形的橫蕩鎖頻。

2)獲取了平臺縱蕩和橫蕩無量綱振幅曲線，發(fā)現(xiàn)鎖頻現(xiàn)象對渦激運動響應幅值影響明顯。

3)物理模型試驗因無法滿足Re相似，高估了黏性力，在橫蕩鎖頻時高估了順流向的運動幅值，低估了橫流向的渦激運動幅值。