胡 彬,水慶象,王大國,彭仁濤

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

特征線算子分裂有限元的圓柱繞流大渦模擬

胡 彬,水慶象,王大國,彭仁濤

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

為研究圓柱繞流流場(chǎng)特性,將大渦模擬與特征線算子分裂有限元相結(jié)合,建立了大渦模擬特征線算子分裂有限元模型,對(duì)單圓柱和串列雙圓柱繞流問題進(jìn)行了數(shù)值模擬,所得結(jié)果與現(xiàn)有研究結(jié)果吻合較好。模擬結(jié)果表明:對(duì)單圓柱繞流,隨著雷諾數(shù)的增大,圓柱近尾流區(qū)上下交替的渦旋逐漸靠近通過圓柱幾何中心的水平線,且渦脫落位置逐漸靠近圓柱。對(duì)Re=1 000的串列雙圓柱繞流,臨界間距在圓柱直徑的2.25～2.5倍之間; 當(dāng)兩圓柱間距小于臨界間距時(shí),上游圓柱后方無明顯渦旋脫落,間隙處壓力較穩(wěn)定;大于臨界間距時(shí),有渦旋脫落,上游圓柱尾流區(qū)上下表面交替出現(xiàn)強(qiáng)負(fù)壓區(qū)。

圓柱繞流;大渦模擬;特征線算子分裂有限元;臨界間距

湍流是流體運(yùn)動(dòng)中的常見現(xiàn)象,隨著計(jì)算機(jī)設(shè)備性能的不斷提高,數(shù)值試驗(yàn)廣泛應(yīng)用于湍流的研究中。目前研究湍流的數(shù)值方法主要有直接數(shù)值模擬[1]、雷諾時(shí)均模擬[2]和大渦模擬(large eddy simulation,LES)[3]3種,其中大渦模擬由于計(jì)算精度較高和計(jì)算量較小而被認(rèn)為是最具潛力的方法之一,其基本思想是將湍流運(yùn)動(dòng)的渦團(tuán)分為大尺度渦和小尺度渦,對(duì)大尺度渦直接數(shù)值求解,小尺度渦對(duì)湍流運(yùn)動(dòng)的影響則通過建立模型求解[4]。

有限元法在數(shù)值模擬方面能表現(xiàn)出良好的幾何邊界和邊界條件適應(yīng)性,并被廣泛用于求解流體動(dòng)力學(xué)問題,但經(jīng)典的Garlerkin有限元法在處理流體對(duì)流問題上容易出現(xiàn)數(shù)值振蕩[5]。為克服該困難,近年來發(fā)展了多種有限元法,目前應(yīng)用較為廣泛的有streamline upwing petrov-Garlerkin (SUPG)法[6]、Taylor-Garlerkin法[7]、特征線-Garlerkin法[8]等。Wang等[9]將Taylor展開引入到特征線-Garlerkin法并結(jié)合算子分裂法的優(yōu)點(diǎn)提出了特征線算子分裂(characteristic-based operator splitting,CBOS)有限元法,該方法將Navier-Stokes方程分裂成擴(kuò)散項(xiàng)和對(duì)流項(xiàng),對(duì)流項(xiàng)借鑒了CBS(characteristic-based split)算法[10]的簡單顯式特征線時(shí)間離散,從而避免了其他有限元法修正權(quán)函數(shù)的困難,王偉等[11]采用該方法對(duì)線性剪切來流的方柱繞流展開了研究。本文建立了大渦模擬CBOS有限元模型,對(duì)單圓柱和串列雙圓柱繞流問題進(jìn)行數(shù)值模擬分析。

1 數(shù)值模型

1.1 大渦模擬控制方程

二維非定常黏性不可壓縮流動(dòng)(忽略能量損失)可由連續(xù)方程和Navier-Stokes方程控制,其無量綱形式為

(1)

(2)

式中:(u1,u2)=(u,v),u為水平方向流速，v為垂向流速;p為壓力;t為時(shí)間;Re為雷諾數(shù);f1、f2分別為水平方向和垂向外力;(x1,x2)=(x,y),x為水平坐標(biāo),y為垂向坐標(biāo)。

采用盒式濾波器對(duì)方程(1)(2)進(jìn)行濾波,得到以下方程[12]:

(3)

(4)

式中:帶上橫杠“-”的變量為過濾后大尺度變量;τij為濾波后產(chǎn)生的亞格子應(yīng)力項(xiàng),采用亞格子渦黏模式表示為

(5)

(6)

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

(7)

為簡潔起見,以下各過濾后大尺度變量均略去上橫杠“-”。

1.2CBOS有限元法

基于CBOS有限元法[9]的基本思想,采用算子分裂法將控制方程(3)和(7)分裂成擴(kuò)散項(xiàng)和對(duì)流項(xiàng):

(8)

(9)

式中:ui,n+θ為擴(kuò)散項(xiàng)(式(8))在n+1時(shí)刻的解,同時(shí)是對(duì)流項(xiàng)(式(9))在n+1時(shí)刻的初值;ui,n+1是對(duì)流項(xiàng)在n+1時(shí)刻的解,同時(shí)也是控制方程(3)和(7)在n+1時(shí)刻的解。

在式(7)中,對(duì)亞格子渦黏系數(shù)υt進(jìn)行線性化處理后,式(7)與式(2)在方程形態(tài)上相同,因此,式(8)與式(9)有限元求解可參考Wang等[9]的有限元求解方法。

1.3 求解步驟

當(dāng)已知n時(shí)刻的值,求解n+1時(shí)刻的值時(shí),求解步驟如下:

a. 以n時(shí)刻的速度場(chǎng)ui,n作為初值,求解方程(6),得到n+1時(shí)刻的亞格子渦黏系數(shù)υt。

b. 以n時(shí)刻的速度場(chǎng)ui,n、壓力場(chǎng)pn作為初值,求解擴(kuò)散方程(8),得到n+1時(shí)刻速度場(chǎng)的過渡值ui,n+θ和壓力值pn+1。

c. 以u(píng)i,n+θ作為初值,求解對(duì)流方程(9),得到n+1時(shí)刻速度場(chǎng)ui,n+1,完成n+1時(shí)刻的求解。

d. 轉(zhuǎn)到下一時(shí)刻,重復(fù)步驟a、b、c。

2 單圓柱繞流模擬

2.1 計(jì)算區(qū)域、邊界條件和網(wǎng)格劃分

2.2 特征參數(shù)

(12)

式中:f0為渦的自然脫落頻率;υ為運(yùn)動(dòng)黏性系數(shù)。

表1 不同Re時(shí)本文模型計(jì)算結(jié)果與其他文獻(xiàn)數(shù)據(jù)比較

從表1可以看出,本文單圓柱繞流特征參數(shù)計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)研究成果接近。

2.3 卡門渦街流動(dòng)特性

圖1為Re=200、1 000、3 900時(shí)半個(gè)渦脫落周期的壓力與流線形態(tài)?？梢钥闯?不同雷諾數(shù)下,圓柱后方渦旋均產(chǎn)生卡門渦街流動(dòng)特性,渦旋上下交替產(chǎn)生、發(fā)展并脫落。另外,隨著雷諾數(shù)的增大,圓柱近尾流區(qū)上下交替的渦旋逐漸靠近通過圓柱幾何中心的水平線,同時(shí),渦脫落位置逐漸靠近圓柱。

圖1 不同Re下半個(gè)渦脫落周期的壓力與流線形態(tài)

3 串列雙圓柱繞流(Re=1 000)模擬

3.1 計(jì)算區(qū)域和網(wǎng)格劃分

圖2為串列雙圓柱繞流計(jì)算域:流場(chǎng)入口距上游圓柱中心5D,出口距下游圓柱中心16D,橫向尺寸為16D,兩圓柱中心間距為L。邊界條件設(shè)置與單圓柱繞流完全相同。共計(jì)算了兩圓柱中心間距L=2D、2.25D、2.5D、4D、5.5D5種工況,圖3為L=2D的網(wǎng)格劃分情況。

圖2 雙圓柱繞流計(jì)算域

圖3 模型網(wǎng)格劃分(L=2D)

3.2 流場(chǎng)特征

圖4為串列雙圓柱繞流在不同間距約半個(gè)周期的流線。當(dāng)L/D≤2.25時(shí)(圖4(a)(b)),從上游圓柱上(下)表面通過的部分流線流至下游圓柱時(shí),回流穿過兩圓柱間隙并從下游圓柱下(上)表面通過,形成兩個(gè)回流區(qū),交替覆蓋下游圓柱。當(dāng)L/D=2.5時(shí)(圖4(c)),流線在兩圓柱間隙發(fā)生“U”形轉(zhuǎn)彎。與間距L/D≤2.25不同的是,一部分之前通過上游圓柱上(下)表面的流線,開始穿過間隙,從下游圓柱的上(下)表面通過。當(dāng)L/D≥2.5時(shí)(圖4(c)(d)),上游圓柱后方出現(xiàn)渦旋脫落。

圖4 不同間距下的流線

可見,Re=1 000的串列雙圓柱繞流,其臨界間距在2.25D～2.5D之間。兩圓柱間距小于臨界間距時(shí),上游圓柱后方無明顯渦旋脫落;大于臨界間距時(shí),有渦旋脫落。下游圓柱始終存在渦旋脫落現(xiàn)象。

3.3 圓柱受力分析

圖5為兩圓柱平均阻力系數(shù)和升力系數(shù)幅值隨L/D的變化曲線,可以看出本文模擬結(jié)果與相關(guān)文獻(xiàn)計(jì)算結(jié)果吻合較好。由圖5(a)可知,隨著L/D的增大,上游圓柱平均阻力系數(shù)逐漸靠近單圓柱繞流計(jì)算結(jié)果;下游圓柱平均阻力系數(shù)在L/D=2.25時(shí)為負(fù),在L/D=2.5時(shí)突然增大為正值,之后減小且趨于穩(wěn)定。由圖5(b)可知,在L/D=2.25～2.5時(shí),上、下游圓柱升力系數(shù)幅值均突然增大；隨后上游圓柱升力系數(shù)幅值逐漸趨近單圓柱繞流計(jì)算結(jié)果，下游圓柱升力系數(shù)幅值穩(wěn)定后大于單圓柱繞流計(jì)算結(jié)果。

圖5 平均阻力系數(shù)和升力系數(shù)幅值隨間距變化的規(guī)律

以上分析表明,上、下游圓柱平均阻力系數(shù)及升力系數(shù)幅值在L/D=2.25～2.5時(shí)均產(chǎn)生突變,進(jìn)一步說明的串列雙圓柱繞流臨界間距在2.25D～2.5D之間,與Jester等[18]的結(jié)果2.38D吻合較好。

圖6為兩圓柱中心間距L=2.25D時(shí),約一個(gè)周期內(nèi)t1=18.47、t2=18.59、t3=18.71、t4=18.95這4個(gè)時(shí)刻的壓力云圖?？梢钥闯?間隙處壓力較穩(wěn)定,上游圓柱上下表面壓力接近,故其升力系數(shù)幅值比單圓柱繞流小很多。同時(shí)間隙處壓力低于下游圓柱近尾流區(qū),故下游圓柱平均阻力系數(shù)為負(fù)值。下游圓柱后方渦脫落位置較遠(yuǎn),對(duì)其背流面上下表面的壓力差影響較小,故升力系數(shù)幅值較單圓柱偏小。

圖6 L=2.25D時(shí)不同時(shí)刻的壓力云圖

圖7 L=2.5D時(shí)不同時(shí)刻的壓力云圖

圖7為兩圓柱中心間距L=2.5D時(shí),約一個(gè)周期內(nèi)t1=9.00、t2=9.12、t3=9.18、t4=9.30、t5=9.42、t6=9.48這6個(gè)時(shí)刻的壓力云圖。可以看出,上游圓柱尾流區(qū)上下表面交替出現(xiàn)強(qiáng)負(fù)壓區(qū),表明有渦旋脫落,故上游圓柱平均阻力系數(shù)及升力系數(shù)幅值接近單圓柱繞流時(shí)的平均阻力系數(shù)及升力系數(shù)幅值。下游圓柱迎流面壓力相對(duì)其尾流區(qū)壓力顯著提高,但依然低于上游圓柱迎流面壓力,故下游圓柱平均阻力系數(shù)較L=2.25D時(shí)增大但小于單圓柱繞流時(shí)的平均阻力系數(shù)。從圖7(a)(b)(c)可以看出,下游圓柱迎流面上表面的正壓力逐漸增加,下表面的負(fù)壓力則逐漸減小,這對(duì)下游圓柱上下表面壓力差產(chǎn)生疊加影響,故下游圓柱升力系數(shù)幅值劇增。

4 結(jié) 論

a. 建立了大渦模擬特征線算子分裂有限元模型,對(duì)單圓柱和串列雙圓柱繞流問題進(jìn)行了數(shù)值模擬,結(jié)果與現(xiàn)有研究結(jié)果吻合較好。

b. 在單圓柱繞流問題中,隨著雷諾數(shù)增大,圓柱近尾流區(qū)交替脫落的渦旋逐漸靠近通過圓柱中心的水平線,且渦脫落位置逐漸靠近圓柱。

c.Re=1 000時(shí)串列雙圓柱繞流臨界間距在2.25D～2.5D之間,兩圓柱平均阻力系數(shù)及升力系數(shù)幅值在臨界間距區(qū)間產(chǎn)生突變。當(dāng)兩圓柱間距小于臨界間距時(shí),上游圓柱后方無明顯渦旋脫落,間隙處壓力較穩(wěn)定;大于臨界間距時(shí),有渦旋脫落,上游圓柱尾流區(qū)上下表面交替出現(xiàn)強(qiáng)負(fù)壓區(qū)。

[1] 朱海濤,單鵬.帶熱傳導(dǎo)效應(yīng)平面渦輪葉珊流動(dòng)直接數(shù)值模擬[J].力學(xué)學(xué)報(bào),2013,45(5):672-680.(ZHU Haitao,SHAN Peng.Direct numerical simulation of turbing cascade flow with heat transfer effects[J].Chinese Journal of Theoretical and Applied Mechanics,2013,45(5):672-680.(in Chinese))

[2] 楊金龍,賴煥新.雷諾時(shí)均模擬的兩種轉(zhuǎn)捩模型比較[J].華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,40(5):524-532.(YANG Jinlong,LAI Huanxin.Numerical investigation of transitional flow by applying two transitional for Reynolds average Navier-Stokes simulation[J].Journal of East China University of Science and Technology(Natural Science Edition),2014,40(5):524-532.(in Chinese))

[3] 胡政敏,肖志堅(jiān),陶鑄,等.基于LES的亞臨界和超臨界圓球繞流數(shù)值模擬[J].水動(dòng)力學(xué)研究與進(jìn)展：A輯, 2016,31(4):496-502.(HU Zhengmin,XIAO Zhijian,TAO Zhu,et al.LES investigations of flow over a sphere in the subcritical and supercritical regimes[J].Chinese Journal of Hydrodynamics,2016,31(4):496-502.(in Chinese))

[4] 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.

[5] 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.

[6] WANG D G,SHUI Q X.SUPG finite element method based on penalty function for lid-driven cavity flow up toRe=27 500[J].Acta Mechanica Sinica,2016,32(1):54-63.

[7] 周桂云,李同春,錢七虎.水庫滑坡涌浪傳播有限元數(shù)值模擬[J].巖土力學(xué),2013,34(4):1197-1201.(ZHOU Guiyun,LI Tongchun,QIAN Qihu.Finite element numerical simulation of water waves due to reservoir landslides[J].2013,34(4):1197-1201.(in Chinese))

[8] ZIENKIEWICZ O C,MORGAN K,SAI B V K,et al.A general algorithm for compressible and incompressible flow:part Ⅱ.tests on the explicit form[J].International Journal for Numerical Methods in Fluids,1995,20(8/9):887-913.

[9] 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.

[10] MALAN A G,LEWIS R W.An artificial compressibility CBS method for modeling heat transfer and fluid flow in heterogeneous porous materials[J].International Journal for Numerical Methods in Engineering,2011,87(1/2/3/45):412-423.

[11] 王偉,水慶象,王大國,等.線性剪切來流對(duì)方柱繞流特性影響的數(shù)值分析[J].水利水電科技進(jìn)展,2017,37(1):17-21.(WANG Wei,SHUI Qingxiang,WANG Daguo,et al.Numerical analysis of influences of linear shear flow on characteristics of flow past a square cylinder[J].Advances in Science and Technology of Water Resources,2017,37(1):17-21.(in Chinese))

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

[13] SMAGORINSKY J.General circulation experiments with the primitive equations:I.the basic experiment[J].Monthly Weather Review,1963,91(3):99-164.

[14] XU S,WANG Z J.An immersed interface method for simulating the interaction of a fluid with moving boundaries[J].Journal of Computational Physics,2006,216(2):454-493.

[15] LIU C,ZHENG X,SUNG C H.Preconditioned multigrid methods for unsteady incompressible flows[J].Journal of Computational Physics,1998,139(1):35-57.

[16] HARICHANDAN A B,ROY A.Numerical investigation of low Reynolds number flow past two and three circular cylinder using unstructured grid CFR scheme [J].International Journal of Heat and Fluid Flow,2010,31:154-174.

[17] 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.

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

[19] 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.

[20] KRAVCHENKO A G,MOIN P.Numerical studies of flow over a circular cylinder atRe=3 900[J].Physics of Fluids,2000,12(2):403-417.

[21] CHEN T Q,ZHANG Q H,CHENG L.Performance investigation of 2D lattice Boltzmann simulation of forces on a circular cylinder [J].Transaction of Tianjin University,2010,16(6):417-423.

[22] BREUER M.Numerical and modeling influences on large eddy simulations for the flow past a circular cylinder[J].International Journal of Heat and Fluid Flow,1998,19(98):512-521.

Large eddy simulation of flow around circular cylinder combined with characteristic-based operator-splitting finite element method

HU Bin, SHUI Qingxiang, WANG Daguo, PENG Rentao

(SchoolofEnvironmentandResources,SouthwestUniversityofScienceandTechnology,Mianyang621010,China)

In order to investigate the characteristics of flow fields around circular cylinders, a numerical model was developed by combing large eddy simulation (LES) with the characteristic-based operator-splitting finite element method. The model was used to simulate flow fields around a circular cylinder and two circular cylinders in tandem arrangement, and the obtained results agree with the existing results. The simulation results show that, for the flow around a circular cylinder, with the increase of the Reynolds number, the vortexes occurring at the top and bottom surfaces alternatively in the near-wake region gradually approach each other and pass along the horizontal line through the geometric center of the cylinder, and the location of vortex shedding gradually approach the cylinder; for the two circular cylinders in tandem arrangement, withRe=1 000, the critical space between the two circular cylinders is 2.25 to 2.50 times their diameters; when the space between the two circular cylinders is smaller than the critical space, vortex shedding does not occur at the rear of the upstream cylinder, and the pressure in the space between the two circular cylinders is stable; when the space between the two circular cylinders is larger than the critical space, vortex shedding occurs, and strong negative pressure regions occur at the top and bottom surfaces alternatively in the near-wake region of the upstream cylinder.

flow around circular cylinder; large eddy simulation; characteristic-based operator-splitting finite element method; critical space

四川省教育廳科研項(xiàng)目(16CZ0013);綿陽市科技計(jì)劃項(xiàng)目(14S-02-6);西南科技大學(xué)研究生創(chuàng)新基金(14YCXJJ0039)

胡彬(1989—),女,碩士研究生,主要從事計(jì)算流體力學(xué)研究。E-mail:hu_bin180@163.com

王大國(1975—),男,教授,博士,主要從事計(jì)算流體力學(xué)研究。E-mail:dan_wangguo@163.com

10.3880/j.issn.1006-7647.2017.02.005

TV131.2;O357

：A

：1006-7647(2017)02-0027-06