基于時(shí)間步長(zhǎng)自適應(yīng)技術(shù)的海上浮基風(fēng)電平臺(tái)繞流數(shù)值模擬分析

王開(kāi)加，程建生，段金輝，王景全，王 亮

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

0 引 言

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

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

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

1.1 概念化設(shè)計(jì)

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

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

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

1.2 參數(shù)的定義

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

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

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

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

1.3 控制方程

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

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

2 數(shù)值計(jì)算

2.1 邊界條件

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

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

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

上下邊界：采用對(duì)稱邊界條件；

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

2.2 離散格式及求解

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

2.3 幾何模型

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

圖3 計(jì)算模型

圖4 網(wǎng)格劃分

3 網(wǎng)格和計(jì)算時(shí)間步長(zhǎng)

3.1 網(wǎng)格劃分

3.2 自適應(yīng)時(shí)間步長(zhǎng)

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

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

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

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

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

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

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

4.1 模擬工況

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

4.2 結(jié)果分析

4.2.1 時(shí)間步長(zhǎng)的影響

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

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

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

圖5 Re=100，α=1.5情況下的時(shí)間步長(zhǎng)時(shí)程曲線

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

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

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

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

圖7 圓柱邊界層附近矢量場(chǎng)

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

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

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

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

5 結(jié) 語(yǔ)

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

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

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

[2] 曹豐產(chǎn)，項(xiàng)海帆. 圓柱非定常繞流及渦致振動(dòng)的數(shù)值計(jì)算[J]. 水動(dòng)力學(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é)報(bào)，1999, 31(1): 100-105.