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

（解放軍理工大學(xué)野戰(zhàn)工程學(xué)院，南京 210007）

0 引 言

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

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

1 數(shù)學(xué)模型

1.1 概念化設(shè)計

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

圖1 海上浮基風(fēng)電平臺概念化設(shè)計

圖2 海上浮基風(fēng)電平臺模型

1.2 參數(shù)的定義

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

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

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

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

1.3 控制方程

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

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

2 數(shù)值計算

2.1 邊界條件

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

進(jìn)口處：給定速度進(jìn)口，u=U,v=0；

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

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

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

2.2 離散格式及求解

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

2.3 幾何模型

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

圖3 計算模型

圖4 網(wǎng)格劃分

3 網(wǎng)格和計算時間步長

3.1 網(wǎng)格劃分

3.2 自適應(yīng)時間步長

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

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

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

為了保證流場在推進(jìn)求解中不發(fā)生“扭曲”現(xiàn)象，選擇在所有網(wǎng)格點（總數(shù)設(shè)為M）上計算(i=1,2,…,M)，從中選擇最小的一個，再附一調(diào)節(jié)系數(shù)α，即?。?/p>

式中：α——自適應(yīng)時間步長適應(yīng)系數(shù)。

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

4 計算結(jié)果及分析

4.1 模擬工況

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

4.2 結(jié)果分析

4.2.1 時間步長的影響

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

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

表1 升阻力系數(shù)及斯特魯哈數(shù)對比表

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

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

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

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

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

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

4.2.2 網(wǎng)格設(shè)計的影響

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

表2 粗糙網(wǎng)格和精細(xì)網(wǎng)格對數(shù)值模擬結(jié)果的影響

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

5 結(jié) 語

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

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

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

[2] 曹豐產(chǎn)，項海帆. 圓柱非定常繞流及渦致振動的數(shù)值計算[J]. 水動力學(xué)研究與進(jìn)展（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] 蘇銘德，康欽軍. 亞臨界雷諾數(shù)下圓柱繞流的大渦模擬[J]. 力學(xué)學(xué)報，1999, 31(1): 100-105.