基于時間步長自適應技術的海上浮基風電平臺繞流數值模擬分析

王開加，程建生，段金輝，王景全，王 亮

（解放軍理工大學野戰工程學院，南京 210007）

0 引 言

隨著能源問題的日趨嚴峻，開發海上風能并向遠海、深海發展是必然趨勢，漂浮式結構是深海基礎的主要方案。而解決圓柱繞流問題是分析深遠海漂浮式風電基礎平臺在波浪和海流的作用下的動力特性的關鍵。對于多圓柱的研究一直是流體力學領域中一個重要研究課題。由于繞流場的位置和形狀的復雜性，數值求解成為十分重要的手段[1]。

目前，對圓柱繞流采用了不同的方法。對不同雷諾數情況下的圓柱繞流問題包括二維數值模擬作了研究。而在二維數值模擬研究中，易獲得較好的阻力系數DC 和斯特魯哈數St，但升力系數LC的數值計算結果與真實值相差較大，這是由于流場網格（特別是邊界層附近區域的流場網格）不夠細密及時間步長的選取不合理造成的。在流固耦合的數值模擬中，網格設計及時間步長是影響數值模擬精度的重要因素。本文在應用自適應時間步長理論及小雷諾數（Re=100）情況下，采用有限體積法對圓柱繞流進行了數值模擬，并對流場特性作了分析。

1 數學模型

1.1 概念化設計

以REpower公司研制開發的5MW的水平軸風力發電機為設計依據，自主設計了一種風電浮式基礎平臺見圖1，結構主體采用半潛式浮式平臺，為了降低整個系統的重心，底部采用懸掛重物的辦法。概念設計的風電浮基系統由3個部分組成：風電機組和塔筒組成的風機系統，甲板、立柱和雙浮體組成的半潛式平臺和包括6根連接桿件在內的壓載物部分。確定了初步設計尺寸參數，并建立了幾何模型，見圖2。

圖1 海上浮基風電平臺概念化設計

圖2 海上浮基風電平臺模型

1.2 參數的定義

在計算流體力學中，一般用無量綱的升力系數LC和阻力系數DC 來衡量某物體的受力情況。升力系數和阻力系數定義如下：

式中：DF和LF——分別為作用于單位長度圓柱上的阻力和升力；U——入流速度。

斯特魯哈數St是描述圓柱繞流的一個非常重要的無量綱量，表征旋渦脫離情況。其定義為：

式中：T——圓柱單側旋渦脫落周期；D——圓柱直徑；U——未受干擾的自由來流速度。

1.3 控制方程

對于不可壓縮黏性流體，在笛卡爾坐標系下，其運動的數值解受控于N-S方程，平面流動的連續性方程和動量方程為：

式中：u、v——分別為x，y方向的速度分量，p——壓力，υ——流體的運動黏性系數，ρ——流體密度。

2 數值計算

2.1 邊界條件

流動方向與垂直流動方向上的均勻速度分別定義為u和v。

進口處：給定速度進口，u=U,v=0；

出口處：在計算域的出口處，設置壓力出口邊界條件；

上下邊界：采用對稱邊界條件；

圓柱表面：給定無滑移、無穿透邊界條件，u=U,v=0。

2.2 離散格式及求解

計算雷諾數（Re=100）為小雷諾數，屬層流范圍，故采用Laminar模型。采用基于壓力基的分離式求解器進行求解，計算中采用具有二階隱式時間格式的非定常流動進行計算。壓力項與速度項的耦合項計算采用PISO算法實現，壓力項離散采用具有二階精度的格式離散，動量方程采用二階迎風格式離散。計算中壓力、密度、體積力和動量項的欠松弛因子分別為0.3、1、1和0.7。

2.3 幾何模型

結合本文研究目標，取圓柱直徑 2mD= ，計算區域為5030DD×的矩形區域，如圖3所示。上游尺寸20D，下游尺寸30D，D為物體垂直于來流方向平面上的特征尺寸，對圓柱一般取直徑。

圖3 計算模型

圖4 網格劃分

3 網格和計算時間步長

3.1 網格劃分

3.2 自適應時間步長

在流固耦合的數值模擬中，時間步長是影響數值模擬精度的重要因素。如果時間步長過大，則圓柱和流場的耦合計算不夠緊密，會導致計算結果精度不高，甚至會出現錯誤的結果。如果時間步長太小，則無謂地增加計算量。目前有關文獻資料中對時間步長的選擇方案有兩種：一種是取圓柱自振周期的0.2%[2]，另外一種是時間步長取遠處來流速度經過圓柱直徑長度所耗費時間(對流時間單位)的0.2%[3]。

在Fluent中，當求解器采用PISO時，即使時間步長取得較大，計算仍然會收斂，但結果的準確性無法得到保證。關于時間步長的選取，引入CFL穩定性限制條件，即

式中：cV——網格單元的體積（二維時取為面積）；d——維數（二維計算取為2，三維計算取為3）；——相應網格單元內流速。式(6)表示計算的時間步長tΔ近似等于流速從一個網格點傳播到另一個網格點所需的時間。由于采用的是非均勻網格以及每個網格點處當地的流速不一樣，所以每個網格點處的tΔ是不一樣的。因此，在給定的時刻t和給定的網格點i上，式(6)應寫成：

為了保證流場在推進求解中不發生“扭曲”現象，選擇在所有網格點（總數設為M）上計算(i=1,2,…,M)，從中選擇最小的一個，再附一調節系數α，即取：

式中：α——自適應時間步長適應系數。

采用上述的自適應時間步長方案，通過用戶自定義功能UDF的二次開發接口對Fluent軟件二次開發來實現自適應時間步長方案在圓柱周圍流場的數值模擬中應用。

4 計算結果及分析

4.1 模擬工況

采用同一條件進行繞流計算：入流速度為1m/s，圓柱直徑D=2m，Re=100。對自適應時間步長適應系數α為0.5、1.0、1.5和固定時間步長為0.025、0.050、0.100共6種工況進行了數值模擬。數值模擬結果見圖5～7，同其他數值模擬結果和實驗結果的比較見表1、2。

4.2 結果分析

4.2.1 時間步長的影響

針對上述6種工況進行數值模擬，給出了采用自適應時間步長時α對計算結果的影響，如表1所示，通過與前人研究結果比較知，該計算結果是可靠的并提高了計算精度。從表1中可以看到，隨著α的增大，St數、升、阻力系數都會逐漸減小。當α=1.0時，St數、升、阻力系數相對于α=1.5時的數據變化不大，但實際運算中計算量卻要大50%，因此全面考慮計算時間、效率以及經驗，取α=1.5時自適應時間步長算法為宜。采用固定時間步長時，St數、升、阻力系數隨時間步長dt增大而減小的趨勢比較明顯。如當dt= 0 .100，計算所得的 St數和阻力系數相比于 dt= 0 .025時減小得并不太多，但升力系數的振幅卻下降了約三分之一。

圖5是Re=100，1.5α=情況下的時間步長時程曲線。從圖中可以明顯地看出，1.5α=情況下采用自適應時間步長算法的時間步長均大于0.025，并且計算結果優于后者。

表1 升阻力系數及斯特魯哈數對比表

圖5 Re=100，α=1.5情況下的時間步長時程曲線

圖6給出了采用自適應時間步長(1.5α=)時圓柱升力系數和阻力系數的時程曲線。由于圓柱尾流形成了卡門渦街，當周期性的渦脫處于圓柱上方時圓柱的升力最大，而后渦脫會逐漸向下游運動，經過圓柱中心時升力系數為0，當渦脫運動到圓柱下游時，升力系數達到最小負值，如此往復，故升力系數的均值為0。阻力系數在渦脫達到穩定后的均值為1.342，這與文獻[4]中Re為100的計算結果非常吻合。

由圖6阻力系數與升力系數時程曲線可以看出，升力系數發生1個周期變化的同時，阻力系數就會發生2個周期的變化，即阻力系數的頻率約為升力系數頻率的2倍。這是由于圓柱發生周期性渦脫時，從上下表面脫離的渦會引起阻力改變1次，而這2個渦共同影響升力變化1次。

圖6 Re=100，α=1.5情況下的圓柱阻力系數與升力系數時程曲線

圖7給出了圓柱邊界層附近的矢量場，從中可以看到在邊界層內，矢量沿圓柱外法線方向逐漸增大，準確地模擬出速度矢量在邊界層內的變化規律。

圖7 圓柱邊界層附近矢量場

4.2.2 網格設計的影響

在靜止圓柱繞流的數值模擬中，比較容易獲得較好的阻力系數DC和斯特哈爾數tS，升力系數LC的數值計算結果有時同真實值相比則偏低。原因是因流場網格劃分不夠精細，特別是邊界層網格劃分不夠精細造成的。在邊界層網格的設計中，除了沿圓周有足夠的網格數，而且第一層網格到圓柱壁面的距離也要足夠小。表2給出了一個比較：

表2 粗糙網格和精細網格對數值模擬結果的影響

從表2的比較可以看出，采用精細網格比采用粗糙網格能獲得更好的數值模擬結果，特別是升力系數的幅值有了提高，更加接近精確值。因此要想獲得較好的數值模擬結果，流場計算區域要足夠大，流場網格劃分要足夠精細，邊界層附近區域的流場網格除了沿圓周要有足夠多的網格數，而且第一層網格到圓柱壁面的距離要足夠小。

5 結 語

深海風電具有風資源豐富、無噪音污染、風電設備利用率高、不受水深和海底地質條件限制等優點，是未來海上風電發展的主要方向。風電浮式基礎，是一種專門為深海區域風電開發的新型海洋基礎，但由于深海區域的各種復雜環境條件引起的作用力對基礎結構的安全造成很大的影響，目前國內外對于深海風電浮基的研究尚處于起步階段。采用時間步長自適應技術和精細的邊界層網格處理技術能夠獲得較為精確的數值結果。

在雷諾數較低、圓柱周圍的流動主要呈現二維流動的情況下，采用基于有限體積法的自適應時間步長算法，流場網格劃分較好、特別是邊界層網格劃分較精細，完全可以用Fluent軟件準確地模擬靜止圓柱繞流問題。本文應用Fluent軟件對單圓柱繞流流場進行數值模擬，所采用的自適應時間步長算法在現有的計算設備下可以獲得較好的數值模擬結果，與固定時間步長算法相比提高了數值計算的精度。為后續風電浮式基礎結構的優化設計和將來深海風電浮基安裝、現場作業等具有指導性的意義。

[1] 王 健，李海濤. 計算流體力學方法在船舶領域的實用性研究[J]. 船舶與海洋工程，2012, (4): 6-11.

[2] 曹豐產，項海帆. 圓柱非定常繞流及渦致振動的數值計算[J]. 水動力學研究與進展（A輯），2001, 16(1): 111-118.

[3] Newman D J, Karniadakis G E. A direct numerical simulation study of flow past a freely vibrating cable[J]. Journal of Fluid Mechanics, 1997, 344: 95-136.

[4] Kim J, Kim D, Choi H. An immersed-boundary finite-volume method for simulations of flow in complex geometries [J]. J Comput Phys, 2001,171:132-150.

[5] Braza M, Chassaing P, Ha Mind H. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder [J]. J Fluid Mech, 1986, 165: 79-130.

[6] Lu XY, Dalton C, Yhang JF. Application of large eddy simulation to an oscillating flow past a circular cylinder [J]. Journal of Fluids Engineering, Transactions of the ASME. 1997, 119(3): 519-525.

[7] 蘇銘德，康欽軍. 亞臨界雷諾數下圓柱繞流的大渦模擬[J]. 力學學報，1999, 31(1): 100-105.