， ，

(哈爾濱工業(yè)大學 能源科學與工程學院,哈爾濱 150001)

1 引 言

鈍體繞流現(xiàn)象是流體力學的經(jīng)典問題之一，國內(nèi)外學者[1-4]對單圓柱、串并列雙圓柱及多個圓柱的繞流進行了深入研究，較好地揭示了其尾跡與受力特性的演化規(guī)律和機制。在實際工程中，除了圓形截面外，具有橢圓形或矩形截面特征的柱體結(jié)構(gòu)也廣泛應用于電力冷卻系統(tǒng)、換熱器和流體機械等領(lǐng)域。其中，橢圓形柱體具有流動阻力低、換熱效果好的優(yōu)點，越來越受到人們的重視，因此開展橢圓柱繞流研究具有重要的實際意義。Dennis等[5]分析了低雷諾數(shù)條件下攻角對單橢圓柱流場結(jié)構(gòu)及升阻力特性的影響規(guī)律；Justin等[6]研究了單橢圓柱繞流尾跡模態(tài)隨長短軸比值的變化情況；Peng等[7]開展了幾何參數(shù)及并列橢圓柱間距對流場尾跡旋渦特性的研究；梁才航等[8]針對橢圓柱流動狀態(tài)及傳熱特性開展了相應的研究工作。上述數(shù)值研究均是基于單個及并列橢圓柱繞流問題展開的，然而在海洋石油平臺支柱和冷凝器排管等實際應用中存在串列的布置形式，因此有必要開展串列橢圓柱繞流的研究工作。

近20多年來，格子Boltzmann方法(Lattice Boltzmann Method)在理論方面取得了很大的進展[9,10]，現(xiàn)已拓展到可圧縮流體、多相流體和燃燒學等應用領(lǐng)域[11-13]，成為計算流體動力學的重要數(shù)值方法。與傳統(tǒng)求解宏觀方程不同，LBM是介于流體的微觀分子動力學模型和宏觀連續(xù)模型之間的介觀模型，主要研究假象粒子的動力學，具有演化過程清晰、邊界處理容易和壓力求解簡單等優(yōu)點[14]。然而，該方法尚存在一些缺點,由于格子速度模型的引用，LBM的求解局限于均勻直角網(wǎng)格，時間推進步長與網(wǎng)格尺度存在相關(guān)性，導致計算較高雷諾數(shù)流動問題時，必須采用較密的網(wǎng)格來保證計算程序的穩(wěn)定性，進而大幅降低了其計算效率。

本文對雷諾數(shù)為100的單橢圓柱及串列雙橢圓柱進行了數(shù)值研究，著重分析了單橢圓柱升阻力系數(shù)和斯特勞哈爾數(shù)隨長短軸比值的變化規(guī)律以及串列雙橢圓柱升阻力和流動形態(tài)隨間距的變化規(guī)律。采用Shu等[15]提出的格子Boltzmann通量求解方法(Lattice Boltzmann Flux Solver)，基于有限體積思想，狀態(tài)變量由宏觀方程給出，界面通量用格子玻爾茲曼的介觀模型構(gòu)造，克服了傳統(tǒng)格子玻爾茲曼方法只能在均勻網(wǎng)格下求解的缺點，并改善了傳統(tǒng)流體動力學計算方法必須花費大量時間求解泊松方程而受到壓力場的限制，大大節(jié)約了數(shù)值計算的時間。此外，采用文獻[16,17]提出的強制浸入邊界法處理橢圓柱物面邊界，以精確滿足無滑移條件。

2 數(shù)值方法

本文采用宏觀不可壓N-S控制方程表達式為

(1)

(2)

在低速流動條件下,基于多尺度Chapman-Enskog展開建立格子波爾茲曼方程中密度分布函數(shù)與宏觀方程狀態(tài)變量及通量之間的關(guān)系，采用浸沒邊界-格子波爾茲曼通量求解(Immersed Boun-dary-lattice Boltzmann Flux Solver) 的不可壓控制方程表達式為

(3)

(4)

(5)

fneqα=-τδt(?/?t+eα·)feqα

(6)

針對以上控制方程，空間離散形式為

(7)

在計算域節(jié)點處理方面，IB -LBFS方法中包括代表流場區(qū)域的歐拉點和代表物理邊界的拉格朗日點，如圖1所示，其本質(zhì)是將浸入邊界模化成流場中的力源項，從而使整個流場的計算可以使用簡單的笛卡爾網(wǎng)格,而不用按照物體表面的形狀生成復雜的貼體網(wǎng)格，簡化了網(wǎng)格生成過程。

計算過程主要分為兩個部分。

(1) 不考慮邊界作用力,使用界面通量重構(gòu)(LBFS)求解控制方程:

(8)

(9)

(2) 采用強制浸入邊界法(IBM)修正宏觀變量，進而求得邊界作用力：

(10)

求解控制方程的關(guān)鍵是求解兩個單元界面上的通量，本文基于格子波爾茲曼模型在網(wǎng)格界面處重構(gòu)數(shù)值格式(LBFS)，如圖2所示，通量表示為

(11)

(12)

通過泰勒級數(shù)展開，非平衡項可以用平衡項表示為

fneqα(r,t)=-τ[feqα(r,t)-feqα(r-eαδt-ri,t)]

(13)

因此只要求得平衡項feqα(r-eαδt-ri,t)與feqα(r,t)便可求得界面通量。

圖1 計算域二維示意圖

Fig.1 A solid boundary immersed in a two -dimensionalcomputational domain

圖2 界面通量重構(gòu)示意圖

Fig.2 Local flux reconstruction at cell interface

已知相鄰兩個網(wǎng)格單元和界面的物理位置，分別定義為ri,ri +1和r，則有

(14)

基于以上方法，可以得到重構(gòu)后界面宏觀變量密度和動量的表達式為

(15)

(16)

再通過LBGK模型求得feqα(r,t)。

在壁面處理方面，應用強制浸入邊界法進行流場速度的修正，使得物面滿足無滑移邊界條件。

定義u*為流場瞬時速度，修正后的歐拉點速度表達式為

un +1=u*+Δu

(17)

歐拉點上修正速度可以通過Dirac delta 函數(shù)插值求解，表達式為

(l=1,…,N) (18)

式中D為連續(xù)kernel函數(shù)，表達式為

(19)

(20)

為了保證無滑移邊界條件，邊界上拉格朗日點的速度等于同一位置流場的速度，則需要通過插值來求解拉格朗日點的流動速度，

(21)

式中l(wèi)=1,…,N;j=1,…,M;h為歐拉點網(wǎng)格尺寸，N和M分別為拉格朗日點與歐拉點的個數(shù)，拉格朗日點上的速度為

(22)

將式(22)寫成矩陣的形式為

AX=B

(23)

(24)

(25)

(26)

3 數(shù)值驗證

(27)

(28)

St=f·D/U∞

(29)

式中Fd和Fl分別為浸入邊界所受到的阻力和升力，fx和fy為邊界點上沿水平和豎直方向上的分力，f為旋渦脫落的無量綱頻率。

由表1數(shù)值計算結(jié)果與其他文獻結(jié)果對比可知，本文計算得到的圓柱繞流升阻力系數(shù)與斯特勞哈爾數(shù)均與文獻結(jié)果吻合良好，充分說明了IB -LBFS求解低雷諾數(shù)不可壓流動繞流問題的精確性。

4 單橢圓柱繞流問題

單橢圓柱繞流的計算域及邊界條件如圖3所示，雷諾數(shù)Re=100，無量綱長度D=Dy=1.0，定義橢圓柱長短軸徑的比值(aspect ratio)為AR=Dx/Dy，保持短軸尺寸不變，通過改變長軸尺寸研究AR=1.0,1.25,1.5,1.75,2.0,2.5,3.0和5.0八種情況下的橢圓柱繞流與受力特性。

從圖4(b)可以看出，最大升力系數(shù)隨AR值的增大而減小，當AR≥3.0時，升力系數(shù)為0，而進一步研究長短軸比值及串列雙橢圓柱間距對流動和受力的影響規(guī)律。圖6為串列雙橢圓柱計算域示意圖，邊界條件與單橢圓柱繞流的邊界條件相同；雷諾數(shù)Re=100，無量綱長度D=Dy=1.0；分別模擬了AR=1.0，1.5，2.0以及L/D=1.0，2.0，3.0，4.0，5.0，7.0，9.0不同間距的流場，其中L代表串列橢圓柱之間的距離。

表1 非定常圓柱繞流升阻力系數(shù)與斯特勞哈爾數(shù)

Tab.1 Drag and lift coefficients，and Strouhal number for unsteady flow past a circular cylinder

ReAuthorsCdClStCalhoun[18]1.35±0.01±0.300.175Russell et al[19]1.38±0.01±0.320.169100Liu et al[20]1.35±0.01±0.340.164Choi et al[21]1.34±0.01±0.320.164Braza et al[22]1.36±0.02±0.340.164Present1.37±0.01±0.340.164Calhoun[18]1.17±0.06±0.670.202Russell et al[19]1.29±0.02±0.500.195200Liu et al[20]1.31±0.05±0.690.192Choi et al[21]1.36±0.05±0.640.191Braza et al[22]1.40±0.05±0.690.192Present1.40±0.05±0.710.192

圖3 橢圓柱繞流計算域

Fig.3 Computational domain of 2D elliptical cylinder

表2表明，斯特勞哈爾數(shù)也隨AR值的增大而減小，當AR≥3.0時，沒有檢測到周期性的尾緣脫落現(xiàn)象，這說明此時橢圓柱繞流為穩(wěn)定流動。

圖5表明隨著AR的增加，橢圓柱繞流時均分離渦與來流垂直方向的尺寸變小，這也說明了其分離角度的減小，而分離渦的長度則先增加后減小。

表2 斯特勞哈爾數(shù)及分離角度隨AR的變化規(guī)律

Tab.2 Strouhal number and separation angle with increasing aspect ratios

圖4 平均阻力系數(shù)及最大升力系數(shù)隨AR的變化規(guī)律

Fig.4 Variation in drag coefficient and maximum lift coefficient withAR

圖5 不同AR下橢圓柱繞流時均流線圖

Fig.5 Time -averaged streamlines with different aspect ratio

5 串列雙橢圓柱繞流問題

圖7給出了串聯(lián)雙橢圓系統(tǒng)中，上下游橢圓柱平均阻力系數(shù)隨間距的變化規(guī)律。串列圓柱阻力系數(shù)與文獻[23,24]的數(shù)值對比再次驗證了本文數(shù)值方法的可靠性，表明IB -LBFS能很好地捕捉層流繞流特性。從上下游橢圓柱的阻力系數(shù)對比可知，在相同AR及L/D的條件下，上游橢圓柱的阻力大于下游橢圓柱阻力，且AR越大，上游橢圓柱阻力系數(shù)越小。

圖7(a)表明，在相同的間距下，上游橢圓柱的阻力系數(shù)隨AR的增加而減小，這與單個橢圓柱阻力的變化規(guī)律相同。在較小的間距下，上游橢圓柱阻力隨著間距的增加呈現(xiàn)小幅度降低的趨勢，這是由于間距的增加使得上游橢圓柱分離的剪切層有較大的空間發(fā)展，分離點后移，引起阻力的降低。隨著間距的繼續(xù)增加，存在一臨界間距，出現(xiàn)阻力躍升現(xiàn)象，當AR=1.0時，該臨界值為L/D=2.0～3.0；當AR=1.5和2.0時，臨界值基本相同，出現(xiàn)在L/D=3.0～4.0之間，此時上游橢圓柱尾跡區(qū)形成脫落渦。當間距大于臨界間距時，橢圓柱阻力變化規(guī)律與圓柱有所不同。對于圓柱，阻力隨間距的變大先變小再變大，逐漸趨于單個圓柱繞流的阻力，這與現(xiàn)有文獻結(jié)論相符；對于橢圓柱，阻力隨間距的變大基本保持不變，再小幅增加，逐漸趨于相同AR值下單個橢圓柱繞流的阻力系數(shù)，說明此時下游橢圓柱對上游橢圓柱受力影響較小。

圖6 串列雙橢圓柱繞流計算域

Fig.6 Computational domain of two tandem elliptic cylinders

圖7(b)表明，在較小的間距下(L/D≤2.0)，AR值越大，下游橢圓柱阻力系數(shù)越大，且AR=1.0 時的阻力系數(shù)為負值，即下游圓柱受到了指向上游的吸力；AR=1.5或2.0時的阻力為正值。當間距達到臨界間距時，下游橢圓柱也出現(xiàn)了阻力躍升。隨著間距進一步增加到L/D=7.0時，下游橢圓柱阻力隨間距的變化規(guī)律略有不同，這是由于不同的間距導致上游橢圓柱的非定常脫落渦結(jié)構(gòu)作用在下游橢圓柱上，使得下游橢圓柱的附著點與分離點位置發(fā)生變化，造成橢圓柱前后壓差以及壁面摩擦力的波動所致。當L/D≥7.0時，阻力系數(shù)隨著間距的增加而緩慢增加，且AR值越大，阻力系數(shù)越大。

為了研究流動形態(tài)的變化，圖8給出了AR=1.0，1.5及2.0情況下，斯特勞哈爾數(shù)隨雙橢圓柱間距的變化規(guī)律以及不同間距下的瞬時渦量圖。

從圖8可以看出，當AR=1.0時，較小間距下St值隨著間距的增加有小幅下降；結(jié)合渦量云圖9可以看出，L/D=1.0和2.0時，上游圓柱剪切層沒有足夠的空間發(fā)展為脫落渦結(jié)構(gòu)而附著在下游圓柱表面上，在下游圓柱后方形成周期性漩渦脫落狀態(tài)，間距的變大使得渦長度增加，因此脫落渦頻率降低；當達到臨界間距時，St突然增大，表明此時上游圓柱尾跡區(qū)域開始形成脫落渦結(jié)構(gòu)；當L/D≥3.0時，St隨間距的增加逐漸趨近于單個圓柱繞流的St值，上游圓柱產(chǎn)生的渦脫渦交替或間歇作用于下游圓柱上，并從下游圓柱表面脫落。由于串列圓柱直徑相同，因此兩圓柱的脫落渦頻率相同。

圖7 橢圓柱阻力系數(shù)隨間距的變化規(guī)律

Fig.7 Dependence of the mean drag coefficient of onthe dimensionless spacingL/D

當AR=1.5和2.0時，在間距1.0≤L/D≤3.0 范圍內(nèi)，流動呈現(xiàn)穩(wěn)定流動狀態(tài)，如圖10和 圖11 所示，在上下游橢圓柱尾跡區(qū)均無脫落渦形成，St值為0，說明此時下游橢圓柱可以有效抑制上游橢圓柱流動分離的現(xiàn)象；當間距達到臨界間距，即L/D=3.0～4.0時，上游橢圓柱尾跡區(qū)域形成的脫落渦結(jié)構(gòu)作用在下游橢圓柱表面，在下游橢圓柱尾跡區(qū)域形成同樣頻率的脫落渦結(jié)構(gòu)，導致St突然增大；隨著間距的繼續(xù)增加，St值逐漸趨于單個橢圓柱繞流的St值。此外，圖10和圖11進一步表明了串列橢圓柱的臨界間距大于串列圓柱。

圖8 斯特勞哈爾數(shù)隨間距的變化規(guī)律

Fig.8 Dependence of the Strouhal number on the dimensionless spacingL/D

圖9AR=1.0時不同間距瞬態(tài)渦量圖

Fig.9 Instantaneous vorticity atAR=1.0

圖10AR=1.5時不同間距瞬態(tài)渦量圖

Fig.10 Instantaneous vorticity atAR=1.5

圖11AR=2.0時不同間距瞬態(tài)渦量圖

Fig.11 Instantaneous vorticity atAR=2.0

6 結(jié) 論

本文應用浸入邊界-格子玻爾茲曼通量求解法對橢圓柱繞流開展了數(shù)值研究。首先通過單圓柱非定常計算驗證了數(shù)值方法的可靠性，進而模擬了低馬赫數(shù)下Re=100時單橢圓及串列雙橢圓柱繞流流動，分析了長短軸比值和間距幾何參數(shù)對其流場及受力特性的影響規(guī)律，結(jié)論如下。

(1) 對于單橢圓柱，平均阻力系數(shù)先隨長短軸比值A(chǔ)R的增加而降低，而后呈現(xiàn)緩慢上升趨勢；升力系數(shù)最大值隨長短軸比值的增大而減小；壓差阻力與摩擦阻力呈相反的變化規(guī)律，當AR<2.5時，壓差阻力起主導作用，而隨著AR值的增加，摩擦阻力起主導作用。隨著AR值的增加，繞流體后方由周期性的旋渦脫落狀態(tài)轉(zhuǎn)變?yōu)榉€(wěn)定對稱流動。

(2) 對于串列雙橢圓柱，上下游橢圓柱阻力隨間距的變化規(guī)律與串列圓柱基本相同,但其臨界間距不同。在較小的間距下，上游橢圓柱阻力隨間距的增加而降低，而下游橢圓柱的阻力隨間距的變大而變大；當間距增大到臨界間距時，上下游橢圓柱阻力均躍升，且串列橢圓柱臨界間距大于串列圓柱，而后隨間距的變大阻力有所波動；當L/D≥7.0時，上下游橢圓柱阻力隨間距增加而增加，逐漸趨近于此長短軸比值下的單橢圓柱阻力。

(3) 在較小的間距下，串列圓柱呈現(xiàn)周期性脫落形式，而串列雙橢圓柱流動狀態(tài)則為穩(wěn)定流動，即下游橢圓柱可以有效地抑制上游橢圓柱流動分離的現(xiàn)象；隨著間距的增加，串列雙橢圓柱與串列圓柱流動狀態(tài)一致，在上下游橢圓柱及圓柱后方均呈現(xiàn)周期性的脫落渦，不斷增加的間距對流動狀態(tài)的影響越來越小，上下游鈍體繞流形式均趨于相同長短軸比值下的單鈍體繞流。

:

[1] Harichandan A B,Roy A.Numerical investigation of low Reynolds number flow past two and three circular cylinders using unstructured grid CFR scheme[J].InternationalJournalofHeatandFluidFlow,2010,31(2):154-171.

[2] Patil D V,Lakshmisha K N.Two -dimensional flow past circular cylinders using finite volume lattice Boltzmann formulation[J].InternationalJournalforNumericalMethodsinFluids,2011,69(6):1149-1164.

[3] Ding H,Shu C,Yeo K S,et al.Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods[J].InternationalJournalforNumericalMethodsinFluids,2006,53(2):305-332.

[4] 張忠宇,姚熊亮,張阿漫.基于間斷有限元方法的并列圓柱層流流動特性[J].物理學報,2016,65(8):236-247.(ZHANG Zhong-yu,YAO Xiong-liang,ZHANG A-man.Numerical simulation of laminar flow past two side -by-side cylinders by discontinuous Galerkin method[J].ActaPhysicaSinica,2016,65(8):236-247.(in Chinese))

[5] Dennis S C R,Young P J S.Steady flow past an elliptic cylinder inclined to the stream[J].JournalofEngineeringMathematics,2003,47(2):101-120.

[6] Leontini J S,Jacono D L,Thompson M C.Stability analysis of the elliptic cylinder wake[J].JournalofFluidMe-chanics,2015,763:302-321.

[7] Peng Y F,Sau A,Hwang R R,et al.Criticality of flow transition behind two side-by-side elliptic cylinders[J].PhysicsofFluids,2012,24(3):034102.

[8] 梁才航,楊永旺,黃斯珉.繞橢圓柱管束的流動與傳熱特性[J].科學技術(shù)與工程.2013,13(13):3592-3597.(LIANG Cai-hang,YANG Yong-wang,HUANG Si-min.Fluid flow and heat transfer across an elliptical cylinder tube bank[J].ScienceTechnologyandEngineering,2013,13(13):3592-3597.(in Chinese))

[9] 郭照立，鄭楚光.格子Boltzmann方法的原理與應用[M].北京：科學出版社，2009.(GUO Zhao -li,ZHENG Chu-guang.TheoryandApplicationsofLatticeBoltzmannMethod[M].Beijing:Science Press,2009.(in Chinese))

[10] 何雅玲，王 勇，李 慶.格子Boltzmann方法的理論及應用[M].北京：科學出版社，2009.(HE Ya-ling,WANG Yong,LI Qing.LatticeBoltzmannMethod:TheoryandApplications[M].Beijing:Science Press,2009.(in Chinese))

[11] Lin C,Xu A,Zhang G,et al.Polar-coordinate lattice Boltzmann modeling of compressible flow[J].PhysicalReviewE,2014,89(1):013307.

[12] Xu A,Lin C,Zhang G,et al.Multiple -relaxation-time lattice Boltzmann kinetic model for combustion[J].PhysicalReviewE,2015,91(4):043306.

[13] Lai H L,Xu A G,Zhang G C,et al.Nonequilibrium thermohydrodynamic effects on the Rayleigh-Taylor instability in compressible flows[J].PhysicalReviewE,2016,94(2):023106.

[14] 聶德明,鄭夢嬌,高 原,等.并列雙橢圓柱繞流的格子Boltzmann虛擬區(qū)域方法的模擬研究[J].計算力學學報,2014,31(4):517-525.(NIE De -ming,ZHENG Meng-jiao,GAO Yuan,et al.Numerical investigation of flow past two elliptical cylinders in side by side arrangement via lattice Boltzmann-fictitious domain method[J].ChineseJournalofComputationalMe-chanics,2014,31(4):517-525.(in Chinese))

[15] Shu C,Wang Y,Teo C J,et al.Development of Lattice Boltzmann flux solver for simulation of incompressible flows[J].AdvancesinAppliedMathematicsandMechanics,2014,6(4):436-460.

[16] Wang Y,Shu C,Teo C J,et al.An immersed boundary-lattice Boltzmann flux solver and its applications to fluid-structure interaction problems[J].JournalofFluidsandStructures,2015,54:440-465.

[17] Wang Y,Shu C.Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications[J].JournalofComputationalPhysics,2009,228(6):1963-1979.

[18] Calhoun D.A Cartesian grid method for solving the two -dimensional streamfunction-vorticity equations in irregular regions[J].JournalofComputationalPhysics,2002,176(2):231-275.

[19] Russell D,Wang Z J.A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow[J].JournalofComputationalPhysics,2003,191(1):177-205.

[20] Liu C,Zheng X,Sung C H.Preconditioned multigrid methods for unsteady incompressible flows[J].JournalofComputationalPhysics,1998,139(1):35-57.

[21] Choi J I,Oberoi R C,Edwards J R,et al.An immersed boundary method for complex incompressible flows[J].JournalofComputationalPhysics,2007,224(2):757-784.

[22] Braza M,Chassaing P,Minh H H.Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder[J].JournalofFluidMechanics,1986,165(1):79-130.

[23] Mussa A,Asinari P,Luo L S.Lattice Boltzmann simulations of 2D laminar flows past two tandem cylinders[J].JournalofComputationalPhysics,2009,228(4):983-999.

[24] Sharman B,Lien F S,Davidson L,et al.Numerical predictions of low Reynolds number flows over two tandem circular cylinders[J].InternationalJournalforNumericalMethodsinFluids,2005,47(5):423-447.