不等直徑串列圓柱繞流大渦模擬

胡 彬，水慶象，王大國

(西南科技大學(xué)環(huán)境與資源學(xué)院，四川綿陽 621010)

不等直徑串列圓柱繞流大渦模擬

胡 彬，水慶象，王大國

(西南科技大學(xué)環(huán)境與資源學(xué)院，四川綿陽 621010)

為研究背負(fù)式海底管線中增設(shè)的小直徑附屬管線對(duì)主管線的水動(dòng)力影響，將大渦模擬中經(jīng)典Smagorinsky亞格子模型與特征線算子分裂有限元法結(jié)合，并引入出口對(duì)流邊界條件，完善了基于特征線算子分裂有限元的大渦模擬方法。通過自編程序數(shù)值模擬Re=1 000的單圓柱繞流，計(jì)算結(jié)果與相關(guān)文獻(xiàn)吻合較好，驗(yàn)證了該算法計(jì)算圓柱繞流的有效性，并分析了Re=1 000時(shí)不同直徑比、間距比情況下的串列雙圓柱繞流，根據(jù)流場的不同渦脫落形態(tài)及兩圓柱平均阻力系數(shù)、升力系數(shù)隨直徑比、間距比變化的規(guī)律得到了不同直徑比條件下的臨界間距范圍。達(dá)到臨界間距后，流場由單一渦脫落狀態(tài)轉(zhuǎn)變?yōu)殡p渦旋脫落狀態(tài)。最后分析了兩圓柱平均阻力系數(shù)及升力系數(shù)在臨界間距后急劇增加的原因，為背負(fù)式海底管線的布局優(yōu)化提供了理論依據(jù)。

大渦模擬； 特征線算子分裂有限元； 不等直徑； 串列； 臨界間距

隨著石油、天然氣的開采逐漸走向深海，在主管線附近鋪設(shè)小直徑附屬管線的背負(fù)式管線[1]在深海油田開發(fā)中顯得日益重要。在背負(fù)式管線中，由于小直徑附屬管線的增設(shè)，主管線附近流場會(huì)發(fā)生較大變化，這種變化對(duì)主管線的受力、尾流結(jié)構(gòu)及渦脫頻率等都將產(chǎn)生影響。因此，對(duì)不等直徑串列雙圓柱繞流的水動(dòng)力研究具有一定的工程意義及應(yīng)用價(jià)值。

研究不等直徑串列圓柱繞流的方法主要有試驗(yàn)方法和數(shù)值模擬方法。Alam & Zhou[2]通過物理試驗(yàn)研究了大圓直徑D=25 mm，小圓直徑d=(0.24～1.00)D，兩圓柱中心間距L=5.5d的雙圓柱繞流。結(jié)果表明：當(dāng)d/D減小時(shí)，小圓柱的斯特勞哈爾數(shù)Sr變小，而大圓柱的Sr數(shù)和平均阻力系數(shù)增大，升力系數(shù)減小。Zhao等[3]用數(shù)值模擬方法研究了直徑比d/D=0.25，L≤1.175D的雙圓柱繞流，結(jié)果表明其尾流結(jié)構(gòu)均為單一渦脫落形態(tài)。Zhao等[4]進(jìn)一步數(shù)值模擬了Re=50 000，d/D=0.5，L≤1.15D的雙圓柱繞流，結(jié)果表明渦脫落形態(tài)仍為單一形態(tài)。Gao等[1]數(shù)值模擬了Re=300時(shí)不等直徑串列雙圓柱繞流，發(fā)現(xiàn)在不同間距條件下出現(xiàn)了單一渦脫落和雙渦旋脫落兩種尾流形態(tài)。于定勇等[5]數(shù)值模擬Re=200時(shí)不同直徑比、間距比情況下的串列雙圓柱繞流，分析了不同直徑比及間距比對(duì)渦脫落形態(tài)、圓柱受力等的影響。

數(shù)值模擬方法研究不等直徑串列圓柱繞流的關(guān)鍵是求解Navier-Stokes(N-S)方程及其湍流模型，其中大渦模擬方法具有較高的計(jì)算精度和較少的計(jì)算量逐漸成為數(shù)值研究復(fù)雜湍流問題的重要方法[6]。由于有限元法具有良好的幾何邊界和邊界條件適應(yīng)性，而被廣泛用于Navier-Stokes(N-S)方程及其湍流模型的求解。但經(jīng)典的Garlerkin有限元法在處理流體對(duì)流占優(yōu)問題上容易出現(xiàn)數(shù)值振蕩[7]。為克服該困難，近年來發(fā)展了多種有限元法，如StreamlineUpwingPetrov-Garlerkin(SUPG)法[8]、Taylor-Garlerkin法[9]、特征線-Garlerkin法[10]。Wang等[11]將Taylor展開引入到特征線-Garlerkin法并結(jié)合算子分裂法的優(yōu)點(diǎn)提出了特征線算子分裂(Characteristic-BasedOperatorSplitting，CBOS)有限元法。該方法將N-S方程分裂成擴(kuò)散項(xiàng)和對(duì)流項(xiàng)，對(duì)流項(xiàng)采用顯式格式，顯式格式具有單個(gè)時(shí)間步計(jì)算量小、程序易于實(shí)現(xiàn)等優(yōu)點(diǎn)。

本文將大渦模擬方法與CBOS有限元法相結(jié)合，數(shù)值模擬了Re=1 000的單圓柱繞流。研究Re=1 000時(shí)不同直徑比、間距比情況下的串列雙圓柱繞流，得到不同直徑比下的臨界間距范圍，分析圓柱受力在臨界間距范圍急劇增大的原因。

1 數(shù)值模型

1.1 大渦模擬控制方程

采用盒式濾波器[12]對(duì)二維非定常無量綱不可壓縮流體的N-S方程組進(jìn)行濾波，根據(jù)湍流動(dòng)能生成與耗散平衡的原則[13]，獲得以下方程組：

?ui/?xi=0

(1)

(2)

式中：帶“-”的變量為過濾后大尺度變量。i,j=1,2；(u1，u2)=(u,v)，u為水平向速度，v為垂向速度；p為壓力；t為時(shí)間；(x1,x2)=(x,y)，x為水平坐標(biāo)，y為垂向坐標(biāo)；雷諾數(shù)Re=Ul/υ(其中，U為特征速度，l為特征長度，υ為運(yùn)動(dòng)黏性系數(shù))。υt稱為亞格子渦黏系數(shù)，Smagorinsky[13]對(duì)其做出如下假設(shè)

(3)

式中：cs為Smagorinsky系數(shù)，一般取0.1～0.2時(shí)可獲得較好的計(jì)算結(jié)果，本文取0.1；Δ為網(wǎng)格過濾尺度。

下面的推導(dǎo)過程為求書寫簡便，略去上橫杠“-”。

1.2CBOS有限元法

在每一個(gè)時(shí)間步內(nèi)，采用算子分裂法將控制方程(1)和(2)分裂成擴(kuò)散項(xiàng)和對(duì)流項(xiàng)：

(4)

(5)

式(4)和(5)應(yīng)用CBOS有限元法[11]，可得如下時(shí)間離散形式：

(6)

(7)

式(7)為對(duì)流項(xiàng)(5)沿特征線顯式時(shí)間離散所得。

圖1 Re=1 000時(shí)單圓柱繞流計(jì)算區(qū)域網(wǎng)格劃分Fig.1 Computation grids of flow past a single cylinder at Re=1 000

如圖1所示，計(jì)算域依左中右及上中下劃分為9個(gè)區(qū)域，由于圓柱附近速度、壓力等梯度大，對(duì)沿柱體徑向及圍繞柱體四周的網(wǎng)格進(jìn)行加密，故每個(gè)區(qū)域的網(wǎng)格尺度各不相同。整個(gè)計(jì)算域劃分為3 238個(gè)9節(jié)點(diǎn)四邊形單元，共有13 212個(gè)節(jié)點(diǎn)數(shù)，其中圓柱表面分布44個(gè)網(wǎng)格，88個(gè)節(jié)點(diǎn)。

表1 計(jì)算的及Sr與其他文獻(xiàn)數(shù)據(jù)的比較

Tab.1ComparisonbetweencalculatedresultswithdatagivenbyotherreferencesatRe=1 000

計(jì)算結(jié)果C-dCALSr本文15712750226HuandKoterayama[14]145—0220JesterandKallinderis[15]1511400250Mittaletal[16]1531370245

2 模型驗(yàn)證

2.1 模型布置

計(jì)算域尺度在主流方向上取為21D，其中圓柱上游分配5D，橫向尺寸為16D，圓柱直徑D=0.1為特征長度；入口處指定沿水平方向的均勻來流U=1為特征速度，垂向速度V=0，Re=1 000。為保證圓柱后方渦旋能夠順利通過出口邊界，出口處為對(duì)流邊界?ui/?t+?ui/?x=0，指定右上角和右下角相對(duì)壓力p=0；側(cè)壁采用可滑移邊界條件；圓柱表面為不可滑移邊界條件。

2.2 特征參數(shù)

圖2 Re=1 000時(shí)單圓柱繞流1個(gè)周期內(nèi)流線Fig.2 Streamline of flow past a single cyliner during a cycle at Re=1 000

3 不等直徑串列雙圓柱繞流分析

3.1 模型布置

圖3(a)給出了串列雙圓柱繞流計(jì)算域：上游圓柱直徑D=0.1為特征長度，d為下游圓柱直徑。流場入口距上游圓柱中心5D，出口距下游圓柱中心16D，橫向尺寸為16D，兩圓柱中心間距為L。取入口處水平方向的均勻來流U=1為特征速度，垂向速度V=0，Re=1 000，指定右上角和右下角相對(duì)壓力p=0；出口處為對(duì)流邊界?ui/?t+?ui/?x=0；側(cè)壁采用可滑移邊界條件；兩圓柱表面為不可滑移邊界條件。

計(jì)算了d/D=0.2，0.4，0.6，0.8，L/D=1.75，1.8，1.9，2.0，2.1，2.2，2.25，2.5共計(jì)32種工況。圖3(b)給出了d/D=0.2，L/D=1.75時(shí)局部網(wǎng)格劃分情況，計(jì)算域依左中右及上中下劃分為12個(gè)區(qū)域，對(duì)柱體徑向及圍繞柱體四周的網(wǎng)格進(jìn)行加密，每個(gè)區(qū)域的網(wǎng)格尺度各不相同。整個(gè)計(jì)算域劃分為4 604個(gè)9節(jié)點(diǎn)四邊形單元，共有18 717個(gè)節(jié)點(diǎn)數(shù)，其中兩圓柱表面均分布40個(gè)網(wǎng)格，80個(gè)節(jié)點(diǎn)。

圖3 不等直徑串列雙圓柱模型布置及網(wǎng)格劃分Fig.3 Computational domain and grid divison when d=0.2D and L=1.75D

3.2 尾流形態(tài)

分別計(jì)算了直徑比d/D為0.2，0.4，0.6，0.8時(shí)在不同間距比下的渦量等值線分布，為節(jié)約版面，只給出d/D=0.2和d/D=0.6時(shí)的渦量等值線分布(見圖4和圖5)，計(jì)算結(jié)果表明，d/D=0.2時(shí)，L/D≤1.8時(shí)上游大直徑圓柱分離的剪切層附著在下游小直徑圓柱表面，渦旋脫落僅發(fā)生在下游圓柱后方，呈現(xiàn)單一渦脫落，渦脫落位置離下游圓柱較遠(yuǎn)；當(dāng)L/D≥1.9時(shí)，上游大圓柱及下游小圓柱均產(chǎn)生渦旋脫落，呈現(xiàn)雙渦旋脫落形態(tài)，渦脫落位置較L/D≤1.8時(shí)更靠近下游圓柱。當(dāng)d/D=0.4，0.6，0.8時(shí)，L/D≤2.1時(shí)上游圓柱后方無渦脫落呈現(xiàn)單一渦脫落形態(tài)；L/D≥2.2時(shí)呈現(xiàn)出雙渦旋脫落流態(tài)。L/D≥2.2時(shí)渦脫落位置離下游圓柱距離大于L/D≤2.1時(shí)的距離。

圖4 d/D=0.2時(shí)不同間距比下的渦量等值線Fig.4 Vorticity contours with different spaces when d/D=0.2

圖5 d/D=0.6時(shí)不同間距比下的渦量等值線Fig.5 Vorticity contours with different spaces when d/D=0.6

3.3 平均阻力系數(shù)及升力系數(shù)幅值

圖6和7分別給出了在不同d/D情況下上、下游圓柱平均阻力系數(shù)及升力系數(shù)隨L/D的變化曲線。由圖可知，d/D=0.2時(shí)，上、下游圓柱平均阻力系數(shù)及升力系數(shù)在L/D=1.8～1.9時(shí)急劇增大；d/D=0.4，0.6，0.8時(shí)，在L/D=2.1～2.2急劇增大。結(jié)合圖4～5流場尾流結(jié)構(gòu)表明：d/D=0.2時(shí)臨界間距范圍為1.8D～1.9D；d/D=0.4，0.6，0.8時(shí)其臨界間距范圍為2.1D～2.2D。

圖6 上游圓柱升、阻力系數(shù)隨L/D的變化曲線Fig.6 Lift and drag coefficients of upper cylinder under different L/D

圖7 下游圓柱升、阻力系數(shù)隨L/D的變化曲線Fig.7 Lift and drag coefficients of lower cylinder under different L/D

3.4 平均阻力系數(shù)及升力系數(shù)急劇變化原因分析

圖8和9分別給出了d/D=0.2， 0.4時(shí)在臨界間距區(qū)間約1個(gè)周期內(nèi)不同時(shí)刻的壓力分布云圖。從圖8(a)和9(a)可以看出：兩圓柱間隙處壓力較穩(wěn)定，上游圓柱上、下表面壓力接近，故其升力系數(shù)比單圓柱繞流小很多。同時(shí)間隙處壓力低于下游圓柱近尾流區(qū)，故下游圓柱平均阻力系數(shù)為負(fù)值。由于未達(dá)到臨界間距，渦脫落位置離下游圓柱較遠(yuǎn)，對(duì)其背流面上、下表面的壓力差影響較小，故升力系數(shù)較單圓柱偏小。

從圖8(b)和9(b)可以看出：達(dá)到臨界間距后，上游圓柱尾流區(qū)上下表面交替出現(xiàn)強(qiáng)負(fù)壓區(qū)，表明有渦旋脫落，故上游圓柱平均阻力系數(shù)及升力系數(shù)值均增大。下游圓柱迎流面壓力相對(duì)其尾流區(qū)壓力顯著提高，故下游圓柱平均阻力系數(shù)均增大。另外，下游圓柱迎流面上表面正壓力逐漸增加(減小)，下表面負(fù)壓力則逐漸減小(增加)，這對(duì)下游圓柱上下表面壓力差產(chǎn)生疊加影響，導(dǎo)致下游圓柱升力系數(shù)振幅劇增。

d/D=0.6，0.8時(shí)上、下游圓柱平均阻力系數(shù)及升力系數(shù)急劇變化的原因與d/D=0.2，0.4時(shí)相近。兩圓柱間距達(dá)到臨界間距后，圓柱周圍壓力變化導(dǎo)致其平均阻力系數(shù)及升力系數(shù)增大。

圖8 d/D=0.2時(shí)臨界間距區(qū)間的壓力分布Fig.8 Pressure distribution at critical spacing range when d/D=0.2

圖9 d/D=0.4時(shí)臨界間距區(qū)間的壓力分布Fig.9 Pressure distribution in critical spacing range when d/D=0.4

4 結(jié) 語

將經(jīng)典Smagorinsky亞格子模型與特征線算子分裂有限元相結(jié)合，通過自編程序模擬了Re=1 000的單圓柱繞流，計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)結(jié)果吻合較好，驗(yàn)證了本文模型計(jì)算圓柱繞流的有效性。研究了Re=1 000，d/D=0.2～0.8，L/D=1.75～2.50的串列雙圓柱繞流，得出如下結(jié)論：

(1)d/D=0.2時(shí)，L/D≤1.8時(shí)渦旋脫落僅發(fā)生在下游圓柱后方，流場呈現(xiàn)單一渦脫落形態(tài)；L/D≥1.9時(shí)，兩圓柱后方均產(chǎn)生渦旋脫落，流場呈現(xiàn)雙渦旋脫落形態(tài)；d/D=0.4，0.6，0.8時(shí)，L/D≤2.1時(shí)流場呈現(xiàn)單一渦脫落形態(tài)，L/D≥2.2時(shí)呈現(xiàn)雙渦旋脫落形態(tài)。

(2)d/D=0.2時(shí)，上、下游圓柱平均阻力系數(shù)及升力系數(shù)在L/D=1.8～1.9時(shí)急劇增大；d/D=0.2，0.6，0.8時(shí)，在L/D=2.1～2.2時(shí)急劇增大。

因此，對(duì)Re=1 000的不等直徑串列雙圓柱繞流：當(dāng)d/D=0.2時(shí)，臨界間距范圍為1.8D～1.9D；當(dāng)d/D=0.4，0.6，0.8時(shí)，其臨界間距范圍為2.1D～2.2D。

[1]GAO Y Y， STEPHONE E， YU D Y， et al. Flow characteristics behind two unequal circular cylinders in tandem arrangement[J]. International Society of Offshore and Polar Engineers， 2010， 1: 1084- 1088.

[2]AlAM M M， ZHOU Y. Strouhal numbers， forces and flow structures around two tandem cylinders of different diameters[J]. Journal of Fluids & Structures， 2008， 24(4): 505- 526.

[3]ZHAO M， CHENG L， TENG B， et al. Numerical simulation of viscous flow past two circular cylinders of different diameters[J]. Applied Ocean Research， 2005， 27(1): 39- 55.

[4]ZHAO M， CHENG L， TENG B， et al. Hydrodynamic forces on dual cylinders of different diameters in steady currents[J]. Journal of Fluids & Structures， 2007， 23(1): 59- 83.

[5]于定勇， 劉洪超， 王昌海. 不等直徑串列雙圓柱體繞流的數(shù)值模擬[J]. 中國海洋大學(xué)學(xué)報(bào)(自然科學(xué)版)， 2012， 42(7/8): 160- 165. (YU Dingyong， LIU Hongchao， WANG Changhai. Numerical simulation of viscous flow past two tandem circular cylinders of different diameters[J]. Periodical of Ocean University of China， 2012， 42(7/8): 160- 165. (in Chinese))

[6]MAHESH K， CONSTANTINESCU G， MOIN P. A numerical method for large-eddy simulation in complex geometries[J]. Journal of Computational Physics， 2004， 197(1): 215- 240.

[7]HACHEM E， RIVAUX B， KLOCZKO T， et al. Stabilized finite element method for incompressible flows with high Reynolds number[J]. Journal of Computational Physics， 2010， 229(23): 8643- 8665.

[8]BROOKS A N， HUGHES T J R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J]. Computer Methods in Applied Mechanics and Engineering， 1990， 32(1/3): 199- 259.

[9]SELMIN V， DONWA J， QUARTAPLLE L. Finite element methods for nonlinear advection[J]. Computer Methods in Applied Mechanics and Engineering， 1985， 52(1/3): 817- 845.

[10]ZIENKIEWICZ O C， MORGAN K， SAI B V K， et al. A general algorithm for compressible and incompressible flow-Part II. Tests on the explicit form[J]. International Journal for Numerical Methods in Fluids， 1995， 20(20): 887- 913.

[11]WANG D G， WANG H J， XIONG J H， et al. Characteristic-based operator-splitting finite element method for Navier-Stokes equations[J]. Science China Technological Sciences， 2011， 54(8): 2157- 2166.

[12]SAGAUT P. Large eddy simulation for incompressible flows[M]. Heidelberg: Springer- Verlag， 2002.

[13]DENHAM， M K， BRIARD P， PATRICK M A. A directionally-sensitive laser anemometer for velocity measurements in highly turbulent flows[J]. Journal of Physics E: Scientific Instruments， 2001， 8(8): 681- 693.

[14]HU C， KOTERAYAMA W. Numerical study on a two-dimensional circular cylinder with a rigid and an elastic splitter plate in uniform flow[J]. International Journal of Offshore and Polar Engineering， 1994， 4(3): 193- 199.

[15]JESTER W， KALLINDERIS Y. Numerical study of incompressible flow about fixed cylinder pairs[J]. Journal of Fluids and Structures， 2003， 17(4): 561- 577.

[16]MITTAL S， KUMAR V， RAGHUVANSHI A. Unsteady incompressible flows past two cylinders in tandem and staggered arrangements[J]. International Journal for Numerical Methods in Fluids， 1997， 25(11): 1315- 1344.

Large eddy simulation of flow past two tandem cylinders with different diameters

HU Bin， SHUI Qingxiang， WANG Daguo

(SchoolofEnvironmentandResources,SouthwestUniversityofScienceandTechnology,Mianyang621010,China)

In order to study the hydrodynamic interaction influences of the small diameter auxiliary pipelines on the main pipelines linked with the piggybacking subsea transport pipeline, a numerical simulation method for the large eddy is proposed on the basis of the characteristic-based operator-splitting finite element method(CBOS), combining the classical Smagorinsky model with the characteristic-based operator-splitting finite element method, and adopting the outlet convective boundary in the numerical simulation. The flow past around a single circular cylinder atRe=1 000 is simulated by the program, and the calculated results well agree with the results given by the other relative literatures, which has validated the efficiency of the calculation method developed by the authors of this paper in simulating the flow past around the circular cylinder. Studies are carried out on the flow past two cylinders atRe=1 000 under the conditions of different diameter ratios and different spaces in a tandem arrangement of two cylinders and the critical spacing ranges with different diameter ratios are obtained, based on the different vortex shedding forms in the flow field and the change characteristics of the mean lift coefficients and the amplitude of the drag coefficients of the large and small cylinders with different diameter ratios and spacing ratios. The wake flow in the flow field indicates that the fluid structure becomes a two-wake shedding mode instead of a single-wake sheddingmode when the gap between the two cylinders is over the critical spacing. In addition, the mean lift coefficients and the amplitude of the drag coefficients of the large and small cylinders both change sharply when the gap is larger than the critical spacing. The analyses of the reasons in the increase of the mean lift coefficients and the amplitude of the drag coefficients of the cylinders are also made when the diameter ratios are equal to 0.2 and 0.4. The research results mentioned above provide a theoretical basis for the layout optimization of the piggybacking propagation subsea pipeline.

large eddy simulation; CBOS finite element method; different diameters; in tandem; critical spacing

10.16198/j.cnki.1009-640X.2017.01.014

2016-01-25

國家自然科學(xué)基金資助項(xiàng)目(41372301)；四川省教育廳科研項(xiàng)目(16CZ0013，15ZB0124)；綿陽市科技計(jì)劃資助項(xiàng)目(14S-02-6)；西南科技大學(xué)研究生創(chuàng)新基金資助項(xiàng)目(14ycx0039)

胡 彬(1989—)，女，四川資陽人，碩士研究生，主要從事計(jì)算流體力學(xué)的研究。E-mail: hu_bin180@163.com 通信作者：王大國(E-mail: dan_wangguo@163.com)

TV131

A

1009-640X(2017)01-0103-08

胡彬， 水慶象， 王大國. 不等直徑串列圓柱繞流大渦模擬[J]. 水利水運(yùn)工程學(xué)報(bào), 2017(1): 103-110. (HU Bin, SHUI Qingxiang, WANG Daguo. Large eddy simulation of flow past two tandem cylinders with different diameters[J]. Hydro-Science and Engineering, 2017(1): 103-110. (in Chinese))