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

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

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

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

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

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

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

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

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

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

1 數值模型

1.1 大渦模擬控制方程

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

(1)

(2)

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

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

(3)

(4)

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

(5)

(6)

式中:cs為Smagorinsky系數,一般取0.1～0.2時可獲得較好的計算結果,本文取0.1;Δ為網格過濾尺度。

(7)

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

1.2CBOS有限元法

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

(8)

(9)

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

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

1.3 求解步驟

當已知n時刻的值,求解n+1時刻的值時,求解步驟如下:

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

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

c. 以ui,n+θ作為初值,求解對流方程(9),得到n+1時刻速度場ui,n+1,完成n+1時刻的求解。

d. 轉到下一時刻,重復步驟a、b、c。

2 單圓柱繞流模擬

2.1 計算區域、邊界條件和網格劃分

2.2 特征參數

(12)

式中:f0為渦的自然脫落頻率;υ為運動黏性系數。

表1 不同Re時本文模型計算結果與其他文獻數據比較

從表1可以看出,本文單圓柱繞流特征參數計算結果與現有文獻研究成果接近。

2.3 卡門渦街流動特性

圖1為Re=200、1 000、3 900時半個渦脫落周期的壓力與流線形態。可以看出:不同雷諾數下,圓柱后方渦旋均產生卡門渦街流動特性,渦旋上下交替產生、發展并脫落。另外,隨著雷諾數的增大,圓柱近尾流區上下交替的渦旋逐漸靠近通過圓柱幾何中心的水平線,同時,渦脫落位置逐漸靠近圓柱。

圖1 不同Re下半個渦脫落周期的壓力與流線形態

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

3.1 計算區域和網格劃分

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

圖2 雙圓柱繞流計算域

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

3.2 流場特征

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

圖4 不同間距下的流線

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

3.3 圓柱受力分析

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

圖5 平均阻力系數和升力系數幅值隨間距變化的規律

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

圖6為兩圓柱中心間距L=2.25D時,約一個周期內t1=18.47、t2=18.59、t3=18.71、t4=18.95這4個時刻的壓力云圖。可以看出,間隙處壓力較穩定,上游圓柱上下表面壓力接近,故其升力系數幅值比單圓柱繞流小很多。同時間隙處壓力低于下游圓柱近尾流區,故下游圓柱平均阻力系數為負值。下游圓柱后方渦脫落位置較遠,對其背流面上下表面的壓力差影響較小,故升力系數幅值較單圓柱偏小。

圖6 L=2.25D時不同時刻的壓力云圖

圖7 L=2.5D時不同時刻的壓力云圖

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

4 結 論

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

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

c.Re=1 000時串列雙圓柱繞流臨界間距在2.25D～2.5D之間,兩圓柱平均阻力系數及升力系數幅值在臨界間距區間產生突變。當兩圓柱間距小于臨界間距時,上游圓柱后方無明顯渦旋脫落,間隙處壓力較穩定;大于臨界間距時,有渦旋脫落,上游圓柱尾流區上下表面交替出現強負壓區。

[1] 朱海濤,單鵬.帶熱傳導效應平面渦輪葉珊流動直接數值模擬[J].力學學報,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] 楊金龍,賴煥新.雷諾時均模擬的兩種轉捩模型比較[J].華東理工大學學報(自然科學版),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] 胡政敏,肖志堅,陶鑄,等.基于LES的亞臨界和超臨界圓球繞流數值模擬[J].水動力學研究與進展：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] 周桂云,李同春,錢七虎.水庫滑坡涌浪傳播有限元數值模擬[J].巖土力學,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] 王偉,水慶象,王大國,等.線性剪切來流對方柱繞流特性影響的數值分析[J].水利水電科技進展,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

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

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

王大國(1975—),男,教授,博士,主要從事計算流體力學研究。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