趙 海 棟， 王 超， 任 冰， 陳 衛(wèi) 東， 王 國 玉

(大連理工大學(xué) 海岸和近海工程國家重點實驗室,遼寧 大連 116024 )

0 引 言

當(dāng)水流經(jīng)過圓柱時，圓柱背流面的兩側(cè)會發(fā)生交替的漩渦脫落，這種周期性的漩渦脫落會使圓柱受到橫流向和順流向流體力作用．該現(xiàn)象廣泛存在于工程實際中，比如海洋石油的輸油管道、鉆井平臺的立管和海底管線等．在實際工程的圓柱結(jié)構(gòu)中，圓柱一般是成組出現(xiàn)，所以雙圓柱繞流問題是國內(nèi)外海洋工程學(xué)者更為關(guān)注的問題．

Meneghini等[1]采用有限單元法對Re=100～200串列和并列布置的兩個圓柱進(jìn)行了數(shù)值模擬．Alam[2]采用有限體積法對Re=200串列雙圓柱進(jìn)行了研究，分析了上下游圓柱漩渦脫落的相位滯后對升力的影響．Zhou等[3]采用有限單元法，使用大渦模擬進(jìn)行紊流計算，對Re=1 000的串列雙圓柱進(jìn)行了數(shù)值模擬．劉松等[4]采用有限體積法對Re=200串列雙圓柱進(jìn)行了研究，分析了不同間距比對圓柱間相互作用和尾流特征的影響．Alam等[5]通過物理模型實驗，得到了亞臨界雷諾數(shù)下串列雙圓柱的受力情況．王龍軍[6]對非等直徑串列雙圓柱的繞流特性進(jìn)行了系統(tǒng)的實驗研究．

然而，目前的研究中少有采用有限差分法(FDM)進(jìn)行圓柱繞流的數(shù)值模擬，原因在于傳統(tǒng)的有限差分法一般采用結(jié)構(gòu)化網(wǎng)格，很難模擬復(fù)雜結(jié)構(gòu)物的邊界．Peskin[7]提出的浸入邊界法(IBM)給這一問題提供了解決思路，該方法通過在結(jié)構(gòu)物附近采用離散的力來模擬結(jié)構(gòu)與流體之間的相互作用，代替流固界面處的物面邊界條件．Fadlun等[8]基于有限差分法和浸入邊界法，提出了一種二階精度、高效的數(shù)值計算方法，該方法可用于解決具有復(fù)雜幾何形狀的物體與三維不可壓縮流體的流固耦合問題．Liu等[9]則進(jìn)一步發(fā)展了這種方法，將浸入邊界法和流體體積法(VOF)結(jié)合獲得虛擬邊界力法(VBF)，再結(jié)合有限差分法，使得該方法在處理包含自由面、復(fù)雜幾何形狀結(jié)構(gòu)物和運動物體的流固耦合問題上具有很大的優(yōu)勢．

此外，不少學(xué)者如Kumar等[10]、Sohankar[11]研究過阻塞率對于圓柱繞流數(shù)值計算結(jié)果的影響，但少有學(xué)者注意到上游圓柱中心與入口的距離對數(shù)值計算結(jié)果的影響，這一參數(shù)足夠大時，才能更好地模擬實際海洋工程中的圓柱繞流問題．

基于以上兩點，展開本文的工作．基于Liu等[9]的工作，應(yīng)用有限差分法和虛擬邊界力法，建立求解均勻流作用下不同間距比的串列雙圓柱繞流的數(shù)學(xué)模型．對數(shù)值計算的控制方程、邊界條件進(jìn)行介紹，同時對計算域大小、網(wǎng)格收斂性和流場三維效應(yīng)進(jìn)行分析；對間距比在2.0～5.0的串列雙圓柱繞流問題進(jìn)行數(shù)值模擬，分析不同間距比下漩渦脫落現(xiàn)象，及上下游圓柱阻力系數(shù)、升力系數(shù)和斯特勞哈爾數(shù)隨間距比變化的趨勢．

1 數(shù)值計算模型

1.1 控制方程

本文采用不可壓縮的連續(xù)方程和黏性的N-S方程作為描述流體運動的基本控制方程：

?·u=0

(1)

(2)

式中：u為流體速度，t為時間，ρ為流體密度，p為流體壓強(qiáng)，ν為運動黏度，g為重力加速度．相較于通用的N-S方程，上式右端增加了一項虛擬邊界力fVBF用于表示結(jié)構(gòu)物對流體的作用力，fVBF的計算見下文．

本文采用經(jīng)典的兩步映射算法對流體的控制方程進(jìn)行數(shù)值離散，兩步映射算法執(zhí)行過程如下：

第一步

(3)

第二步

(4)

(5)

(6)

數(shù)值模型基于有限差分法，采用非均勻交錯網(wǎng)格劃分計算域，其中矢量(如作用力、速度等)定義在網(wǎng)格面心上，標(biāo)量(如壓強(qiáng)、密度等)定義在網(wǎng)格體心上．采用虛擬邊界力法處理結(jié)構(gòu)物邊界，實現(xiàn)結(jié)構(gòu)物和流場之間的數(shù)據(jù)交互．由于傳統(tǒng)的一階迎風(fēng)格式會在計算過程中引入過多的數(shù)值耗散，而中心差分格式常會導(dǎo)致數(shù)值計算的不穩(wěn)定，本文中對流項采用迎風(fēng)格式與中心差分格式相結(jié)合的混合差分格式進(jìn)行離散，可以在減少數(shù)值耗散的同時提高計算的穩(wěn)定性，擴(kuò)散項則采用中心差分格式離散．

1.2 邊界條件

流域邊界條件按物理量分為速度邊界和壓強(qiáng)邊界，如下：

(1)速度入口邊界設(shè)置為均勻流，即u=0.01 m/s，v=0 m/s；出口邊界設(shè)置為開放邊界，假定出口處流體流動已經(jīng)衰退為均勻流動狀態(tài)，速度沿出口法線方向梯度為零，?u/?x=0，?v/?x=0；上下邊界設(shè)置為自由滑移邊界條件，即?u/?y=0，v=0 m/s；流固界面條件設(shè)置為無滑移和不可穿透固體邊界．

(2)壓強(qiáng)入口邊界、出口邊界及上下邊界均設(shè)置為零梯度邊界，沿邊界法線方向梯度為零．

1.3 數(shù)值計算域

為保證數(shù)值計算的精度和效率，對Re=200串列雙圓柱的數(shù)值計算域進(jìn)行討論，并與Alam[2]的計算結(jié)果進(jìn)行對比．選擇圓柱中心間距比L*=4的工況：圓柱繞流數(shù)值模擬采用雷諾相似，上下游圓柱直徑D均為0.1 m，均勻來流速度u0為0.01 m/s，由Re=u0D/ν經(jīng)計算得運動黏度ν=5.0×10-6s/m2．

(1)計算域收斂性驗證

串列雙圓柱的計算域布置見圖1．

圖1 串列雙圓柱的計算域Fig.1 Computational domain for two tandem cylinders

表1 串列雙圓柱的計算域收斂性分析Tab.1 Computational domain convergence analysis for two tandem cylinders

表1中序號1、2所代表的計算域均為40D×30D，此時只改變Lu、Ld的相對取值，計算結(jié)果相差很大，最大為10.2%(Cl1，max)，最小也有4.2%(Cd1，mean)；序號1、3所代表的計算域中Ld=31D，此時只改變Lu的取值(5D、15D)，計算結(jié)果相差很大，最大為11.6%(Cl2,max)，最小也有5.7%(St)；可見Lu的取值對計算結(jié)果影響很大，原因在于來流為均勻流，至上游圓柱前端速度降為零，有較大的速度梯度，需要足夠的網(wǎng)格解析度和空間來捕捉流場演變．3、4號計算域與5號計算域的計算結(jié)果除Cd2，mean(7.5%，7.0%)外，相對差值均在5%以內(nèi)．6號計算域65D×30D與文獻(xiàn)[2]的計算域完全一致，兩者計算結(jié)果相對差值最大值只有2.4%(Cd1，mean)，而與5號計算域的計算結(jié)果相對差值最大值有5.4%(Cd1，mean)．本文后續(xù)計算均選取最大的5號計算域100D×80D，以保證入口處的均勻來流以及減少阻塞率對于流體流速的影響，精度較高．同時本文采用等比格式劃分網(wǎng)格，邊緣區(qū)域網(wǎng)格較大，故加大計算域，網(wǎng)格數(shù)量增加相對較少，從而對計算效率影響較小．

(2)網(wǎng)格收斂性驗證

為提升計算效率，避免資源浪費，同時保證數(shù)值計算精度，進(jìn)行網(wǎng)格收斂性驗證，計算結(jié)果見表2．最小網(wǎng)格1/50D、1/40D、1/30D的計算結(jié)果相差均在1.5%以內(nèi)，所以最小網(wǎng)格選擇1/30D，這樣既可保證計算精度，又能提高計算效率．

表2 串列雙圓柱的網(wǎng)格收斂性分析Tab.2 Grid convergence analysis for two tandem cylinders

(3)二維、三維數(shù)值模擬結(jié)果對比

為確定Re=200串列雙圓柱的適用模型維度，進(jìn)行了三維數(shù)值模擬．圓柱長度設(shè)置為3.2D，網(wǎng)格數(shù)為32，計算結(jié)果為Cd1，mean=1.235，Cl1，max=0.719，Cd2，mean=0.347，Cl2，max=1.512，St=0.170．從圖2可知，Re=200的渦量場沒有明顯的三維現(xiàn)象．對比表2中二維數(shù)值模擬的計算結(jié)果差別極?。示湍壳暗挠嬎阈枨髞碚f，本文使用二維數(shù)值模型進(jìn)行計算．圖2中渦量的計算公式為Ωz=?v/?x-?u/?y，渦量使用均勻來流速度u0和圓柱直徑D進(jìn)行量綱一化．

圖2 串列雙圓柱的三維瞬時渦量場Fig.2 The three-dimensional instantaneous vorticity field of two tandem cylinders

2 數(shù)值計算結(jié)果與分析

使用建立好的數(shù)值模型模擬了Re=200、間距比L*=2.0～5.0的串列雙圓柱的圓柱繞流，將數(shù)值計算結(jié)果與文獻(xiàn)[2]進(jìn)行了比較，結(jié)果詳見表3．

從表3可以看出，本文與文獻(xiàn)[2]在各系數(shù)的計算中趨勢一致，絕對差值均在0.043以內(nèi)．在Cd1，mean計算中，本文與文獻(xiàn)[2]相對差值均在±3%以內(nèi)；在Cd2，mean計算中，本文與文獻(xiàn)[2]相對差值較大，在L*=3.5時達(dá)到了95.56%，但由于基數(shù)較小，絕對差值只有0.043，綜合考慮在可接受范圍內(nèi)；St的相對差值只有在L*=2.0時達(dá)到了6.35%，其余情況相對差值均在±4%以內(nèi)．Wang等[12]也進(jìn)行過L*=2.0的串列雙圓柱的繞流計算，其中Cd1，mean=1.039，Cd2，mean=-0.198，St=0.133，與本文計算結(jié)果更加相符．計算結(jié)果差值產(chǎn)生的原因是本文采用的計算域100D×80D，相較于文獻(xiàn)[2]的計算域65D×30D更大，更能保證入口是均勻來流，以及減少阻塞率對流體流速的影響，這在計算域收斂性驗證中已經(jīng)進(jìn)行過討論，另一方面原因來自于數(shù)值計算方法的差異．

表3 串列雙圓柱在不同間距比下的計算結(jié)果Tab.3 The computational result of two tandem cylinders at different spacing ratios

串列雙圓柱各間距比的瞬時渦量場見圖3．由圖3可知，當(dāng)間距比為2.0、3.0、3.5時，為再附著狀態(tài)(reattachment)，上游圓柱分離的剪切層會重新附著到下游圓柱上，上下游圓柱之間會形成對稱的漩渦．當(dāng)間距比為4.0、4.5、5.0時，為同脫落狀態(tài)(co-shedding)，上下游圓柱均有漩渦脫落．從再附著狀態(tài)到同脫落狀態(tài)轉(zhuǎn)變的間距比稱為臨界間距比Lcr．

為進(jìn)一步研究Re=200串列雙圓柱繞流的臨界間距比Lcr，本文又計算了間距比為3.6、3.7、3.8、3.9的工況，與表3的數(shù)值計算結(jié)果一同繪制在圖4～7中，可見臨界間距比約為3.7．

圖4、5分別展示了上下游圓柱平均阻力系數(shù)隨間距比L*變化的趨勢．當(dāng)L*<3.7時，Cd1，mean

(a)L*=2.0

圖4 上游圓柱在不同間距比下的平均阻力系數(shù)Fig.4 The mean drag coefficients of the upstream cylinder at different spacing ratios

圖5 下游圓柱在不同間距比下的平均阻力系數(shù)Fig.5 The mean drag coefficients of the downstream cylinder at different spacing ratios

隨著L*的增大而減小，Cd2，mean隨著L*的增大而增大，這與剪切層附著點有關(guān)．L*>3.7時，上下游圓柱均有漩渦脫落，Cd1，mean和Cd2，mean均急劇增加后緩慢增加；隨著間距比的增大，上游圓柱阻力系數(shù)受下游圓柱的影響越來越小，趨向于單圓柱狀態(tài)，而下游圓柱處于上游圓柱的尾流中，其阻力系數(shù)雖有所增大，但仍保持較小值．

圖6展示了斯特勞哈爾數(shù)隨間距比變化的趨勢．當(dāng)L*<3.7時，隨著間距比的增大，上游圓柱分離的剪切層需要更多的時間重新附著到下游圓柱，然后從下游圓柱脫落，完成漩渦脫落的循環(huán)，進(jìn)而導(dǎo)致St逐漸減?。?dāng)L*>3.7時，上下游圓柱均有漩渦脫落，St激增至0.175后，隨著間距比的增大，緩慢逼近單圓柱的St．

圖6 串列雙圓柱在不同間距比下的斯特勞哈爾數(shù)Fig.6 The Strouhal number of two tandem cylinders at different spacing ratios

圖7展示了本數(shù)值模型中上下游圓柱的升力系數(shù)最大值隨間距比變化的趨勢，水平線代表單圓柱繞流的升力系數(shù)最大值．在L*<3.7時，上游圓柱沒有漩渦脫落，升力系數(shù)幾乎為零，且隨間距比的增大基本保持不變；隨著間距比的增大，上游圓柱分離的剪切層至下游圓柱時，流速恢復(fù)到更高的水平，導(dǎo)致下游圓柱升力系數(shù)緩慢上升，在L*=3.7時達(dá)到單圓柱升力系數(shù)的55.8%．在L*>3.7時，上下游圓柱均有漩渦脫落，升力系數(shù)均急劇增長；隨著間距比的增大，上游圓柱受下游圓柱的影響越來越弱，升力系數(shù)基本與單圓柱升力系數(shù)持平；下游圓柱在上游圓柱的尾流中，不僅受自身瀉渦的影響，還受上游圓柱瀉渦的影響，升力系數(shù)急劇增長為單圓柱升力系數(shù)的2倍左右，但隨著間距比的增大，上游圓柱的瀉渦至下游圓柱前端會慢慢衰退，逐漸恢復(fù)至均勻流，對下游圓柱的影響越來越小．可以預(yù)見，超過某一間距比后，上下游圓柱均表現(xiàn)得與單圓柱無異．

圖7 上下游圓柱在不同間距比下的升力系數(shù)最大值Fig.7 The maximum lift coefficients of upstream and downstream cylinders at different spacing ratios

3 結(jié) 論

(1)在對計算域收斂性驗證的過程中發(fā)現(xiàn)，上游圓柱中心與入口的距離對數(shù)值模擬的結(jié)果具有較大影響，基于本文數(shù)值結(jié)果，這個距離要在15D以上才能保證數(shù)值計算的精度．

(2)Re=200條件下，三維數(shù)值模型中渦量場不存在明顯的三維特性，上下游圓柱阻力系數(shù)、升力系數(shù)及斯特勞哈爾數(shù)的計算結(jié)果與二維數(shù)值模型基本一致，可以使用二維數(shù)值模型進(jìn)行計算．

(3)Re=200條件下，串列雙圓柱的漩渦脫落狀態(tài)從再附著狀態(tài)到同脫落狀態(tài)轉(zhuǎn)變的臨界間距比為3.7．

(4)低于臨界間距比時，上游圓柱分離的剪切層重新附著到下游圓柱上，上下游圓柱之間會形成對稱的漩渦，漩渦只從下游圓柱脫落．隨著間距比的增大，Cd1，mean和St逐漸減小，Cd2，mean和Cl2，max逐漸增大，Cl1，max幾乎保持不變．

(5)超過臨界間距比時，上下游圓柱均有漩渦脫落，上下游圓柱阻力系數(shù)、升力系數(shù)及斯特勞哈爾數(shù)均有突增．隨著間距比的增大，上游圓柱阻力系數(shù)和升力系數(shù)受下游圓柱的影響越來越小，趨向于單圓柱狀態(tài)；下游圓柱處于上游圓柱的尾流中，其阻力系數(shù)雖有所增大，但仍保持較小值，升力系數(shù)相當(dāng)于單圓柱升力系數(shù)的2倍左右，同時隨著間距比的增大，有下降的趨勢．可以預(yù)見，超過某一間距比后，上下游圓柱均表現(xiàn)得與單圓柱無異．