蘭小杰，趙偉文，萬德成，鄒 麗

(1.上海交通大學(xué) 船海計(jì)算水動(dòng)力學(xué)研究中心(CMHL) 船舶海洋與建筑工程學(xué)院，上海 200240;2.大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024)

柱體繞流問題是流體力學(xué)領(lǐng)域的經(jīng)典問題和研究熱點(diǎn)之一。當(dāng)流體流過柱體時(shí)，會(huì)發(fā)生分離流、渦流等很多復(fù)雜的流動(dòng)現(xiàn)象，準(zhǔn)確地模擬這些現(xiàn)象一直是計(jì)算流體力學(xué)的目標(biāo)之一。同時(shí)柱體繞流也是現(xiàn)代工程中經(jīng)常遇到的問題，風(fēng)吹過高層建筑物，水流流過橋墩，海流流過海洋平臺(tái)等都是典型的柱體繞流。當(dāng)流體流過這些結(jié)構(gòu)物時(shí)，結(jié)構(gòu)物會(huì)受到很大的作用力，對(duì)結(jié)構(gòu)的強(qiáng)度提出了更高的要求。

柱體繞流經(jīng)常被作為標(biāo)準(zhǔn)驗(yàn)證算例，對(duì)其的研究已經(jīng)進(jìn)行得相當(dāng)廣泛。Perry等[1]、Davis等[2]、Okajima等[3]許多學(xué)者用試驗(yàn)的方式對(duì)方柱繞流的流場(chǎng)及方柱受力進(jìn)行測(cè)量，還有大量學(xué)者用基于網(wǎng)格的計(jì)算流體力學(xué)方法對(duì)方柱或圓柱繞流進(jìn)行了相關(guān)研究[4-8]。但是利用無網(wǎng)格粒子法對(duì)柱體繞流進(jìn)行模擬的研究還不算很多。Yildiz等[9]提出了一種用于光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics method，簡(jiǎn)稱SPH)方法的改進(jìn)固體邊界處理方法，多邊界切線法(multiple boundary tangent method，簡(jiǎn)稱MBT)解決了粒子法不能正確處理復(fù)雜流域中彎曲邊界的問題，第一次用SPH方法模擬了雷諾數(shù)為50的圓柱繞流。Shadloo等[10]用不可壓縮光滑粒子法(ISPH)模擬雷諾數(shù)為300的方柱繞流，說明了ISPH方法能夠自然地捕獲鈍體繞流的流分離、渦旋脫落等復(fù)雜的物理現(xiàn)象。Marrone等[11]使用弱可壓縮的SPH方法模擬了低雷諾數(shù)時(shí)鈍體周圍黏性流的演變。Koshizuka等[12]用MPS方法對(duì)雷諾數(shù)為100和200的柱體繞流進(jìn)行了模擬，成功地模擬出了存在于圓柱體和方柱體后面的卡門渦街。國(guó)內(nèi)孫鵬楠[13]對(duì)傳統(tǒng)光滑粒子法進(jìn)行改進(jìn)，并用改進(jìn)后的方法對(duì)圓柱繞流進(jìn)行模擬，對(duì)模擬得到的流場(chǎng)進(jìn)行細(xì)致的分析并與試驗(yàn)結(jié)果進(jìn)行對(duì)比，證明光滑粒子法模擬圓柱繞流是有效的。但相較于圓柱，方柱由于尖銳直角的存在，對(duì)其進(jìn)行模擬更具挑戰(zhàn)性。

MPS方法是在SPH方法的基礎(chǔ)上發(fā)展起來的，最先由Koshizuka和Oka[14]于1996年提出，主要用于求解不可壓縮流動(dòng)問題。MPS 繼承了SPH 的無網(wǎng)格思想，用拉格朗日粒子攜帶空間流場(chǎng)的信息，粒子具有質(zhì)量、動(dòng)量、能量等物理量，粒子間的影響通過核函數(shù)來實(shí)現(xiàn)。近年來，MPS方法得到了長(zhǎng)足的發(fā)展，除了在液艙晃蕩[15]等粒子法傳統(tǒng)優(yōu)勢(shì)領(lǐng)域的應(yīng)用外，研究人員還把MPS方法的應(yīng)用擴(kuò)展到了多相流動(dòng)[16-17]、聲學(xué)[18]等領(lǐng)域。采用基于MPS方法自主開發(fā)的無網(wǎng)格法求解器MLParticle-SJTU，對(duì)雷諾數(shù)分別為40、200及1 000的方柱繞流進(jìn)行仿真模擬分析，得到3種雷諾數(shù)下方柱繞流流場(chǎng)及方柱受力，并與文獻(xiàn)中的結(jié)果進(jìn)行對(duì)比驗(yàn)證。

1 數(shù)值方法

1.1 控制方程

MPS方法的控制方程包括連續(xù)性方程和N-S方程，對(duì)于不可壓縮流體，分別可以寫成式(1)、式(2)形式：

(1)

(2)

式中：P、V分別為壓力和速度，ν、ρ分別表示流體的運(yùn)動(dòng)黏性系數(shù)和密度，f表示質(zhì)量力。MPS方法的控制方程為拉格朗日形式，因此不存在對(duì)流項(xiàng)。

1.2 核函數(shù)

與SPH方法不同，MPS中核函數(shù)只是起著權(quán)函數(shù)的作用，求解過程不需要使用核函數(shù)的導(dǎo)函數(shù)，所以MPS方法的核函數(shù)只要求連續(xù)不要求光滑。采用的核函數(shù)選用Zhang等[19]推薦的核函數(shù)：

(3)

其中，l=|ri-rj|，表示粒子間距，le表示粒子的影響半徑，一般根據(jù)需要選取合適的值即可。

1.3 不可壓縮條件

采用被Lee等[20]改寫過的混合源項(xiàng)法，表達(dá)式為：

(4)

(5)

其中，r是粒子的坐標(biāo)，i、j為粒子的編號(hào)。

1.4 梯度模型

采用的梯度模型表達(dá)式為：

(6)

1.5 散度模型

在MPS方法的連續(xù)性方程中存在散度項(xiàng)，需要用核函數(shù)對(duì)其進(jìn)行離散。采用的散度模型表達(dá)式為：

(7)

1.6 Laplacian模型

在文中MPS方法里使用的是Koshizuka[14]所給出的Laplacian模型，模型表達(dá)式為：

(8)

其中，λ是修正系數(shù)，引入λ是為了修正數(shù)值結(jié)果，使其與擴(kuò)散方程的解析結(jié)果相一致，其表達(dá)式為：

(9)

2 數(shù)值模擬

2.1 計(jì)算模型

研究中，計(jì)算域示意如圖1所示。設(shè)定方柱邊長(zhǎng)D為0.1 m，計(jì)算域大小設(shè)置為25D×12D。為消除入口邊界對(duì)柱體繞流模擬的影響，計(jì)算域流體入口距柱體中心7D，計(jì)算域出口與柱體中心的距離為18D，從而確保尾流得到充分發(fā)展。為減少壁面邊界對(duì)柱體周圍流場(chǎng)的影響，兩側(cè)壁面距柱體中心的距離均設(shè)置為6D。計(jì)算坐標(biāo)系原點(diǎn)在計(jì)算域左下角，x軸正方向?yàn)閬砹魉俣确较颉ＡＷ娱g距設(shè)置為Δx/D=25，計(jì)算雷諾數(shù)分別為Re=40、200及1 000。

圖1 計(jì)算域示意Fig.1 Computational domain of flow around square

2.2 邊界條件

2.2.1 入口邊界條件

圖2為流入流出邊界的示意。圖2(a)為初始狀態(tài)，流入邊界上的壁面粒子以指定的流入速度向右移動(dòng)，流體顆粒被左側(cè)壁面推動(dòng)，也會(huì)以指定的速度向右流動(dòng)，如圖2(b)。當(dāng)壁面向右移動(dòng)距離達(dá)到指定值時(shí)，壁面退回初始位置，原推板位置填入虛粒子(圖2(c))，賦予原推板位置處的虛粒子速度和壓力等物理量，將虛粒子轉(zhuǎn)變成為流體粒子，計(jì)算繼續(xù)進(jìn)行(圖2(d))。推板不斷往復(fù)運(yùn)動(dòng)，推動(dòng)流場(chǎng)中所有流體粒子以指定速度向右流動(dòng)，完成流動(dòng)的模擬。

圖2 流入流出邊界示意Fig.2 Diagrammatic drawing of inlet and outlet

2.2.2 出口邊界條件

如圖2所示，當(dāng)流體粒子流出計(jì)算域時(shí)，把流體粒子轉(zhuǎn)換成虛粒子，不影響其它流體粒子。虛粒子的質(zhì)量、密度等物理量都是0，等待替補(bǔ)到入口邊界。虛粒子填入由于推板后退產(chǎn)生的空隙中后，會(huì)被賦予相應(yīng)的物理性質(zhì)，變成新的流體粒子，重新進(jìn)入計(jì)算域。

2.2.3 壁面邊界條件

在MPS方法中，壁面邊界由多層粒子組成，如圖2中所示。內(nèi)層與流體顆粒相鄰的邊界粒子稱為第一類邊界粒子，第一類邊界粒子參與壓力Poisson 方程的求解。外層由第二類邊界粒子組成，第二類邊界粒子的壓力通過周圍的流體粒子和第一類邊界粒子外插出來。柱體壁面邊界粒子除了速度始終設(shè)為0以外，其它屬性如作用域半徑、密度和質(zhì)量等與流體粒子完全相同，柱體壁面為無滑移壁面。上下邊界壁面給滑移壁面條件，可以更好地模擬均勻來流。

2.3 來流條件

文中來流條件設(shè)定為均勻來流，通過模擬管道內(nèi)均勻流對(duì)均勻來流條件進(jìn)行驗(yàn)證。選取柱體繞流雷諾數(shù)Re=40的工況，設(shè)定來流速度為0.000 4 m/s，模擬沒有柱體時(shí)的管道流動(dòng)。圖3為模擬得到的整體速度云圖，管道整體速度都在0.000 4 m/s附近。圖4為流動(dòng)充分發(fā)展后x/D=7、9、12、17、22處的x方向時(shí)均速度，其中x/D=7處是放置方柱的位置。從圖4中可以看見，不同位置下的x方向速度剖面都呈直線且在0.000 4 m/s 處，誤差小于0.1%，精確模擬均勻管道流。

圖4 不同流向位置處x方向速度Fig.4 The velocity in tne x direction at different positions

3 結(jié)果分析

3.1 水動(dòng)力系數(shù)

在方柱繞流問題中，方柱主要受到的水動(dòng)力包括阻力和升力。一般為了方便對(duì)比研究，需要對(duì)方柱受力進(jìn)行無量綱化處理，得到時(shí)均阻力系數(shù)Cd、升力系數(shù)Cl和斯特勞哈爾數(shù)St，分別為：

(10)

式中：Fd和Fl分別為圓柱受到的阻力和升力；ρ為流體密度；U為流速；fs為斯特勞哈爾渦泄頻率。表1為文中模擬計(jì)算得到的方柱繞流水動(dòng)力系數(shù)結(jié)果。

表1 阻力系數(shù)和斯特勞哈爾數(shù)Tab.1 Drag coefficient and Strouhal number

以Sen等[21]的數(shù)值計(jì)算結(jié)果為參考，比較雷諾數(shù)Re=40時(shí)方柱的阻力系數(shù)Cd。文中計(jì)算得到Re=40時(shí)方柱的阻力系數(shù)為1.724，略大于Sen等[21]的計(jì)算結(jié)果1.67，誤差為3.23%。平均升力系數(shù)為-0.009 2，接近于0，與試驗(yàn)結(jié)果及理論結(jié)果相符合。圖5和圖6分別是Re=40時(shí)升力系數(shù)時(shí)歷曲線和升力功率譜密度，此時(shí)的升力系數(shù)曲線沒有呈現(xiàn)出周期性，其中升力功率譜密度縱坐標(biāo)|P1(f)|表示升力在各頻率下的功率分量大小。

圖5 Re=40，阻力系數(shù)和升力系數(shù)時(shí)歷曲線Fig.5 Time history of the drag and lift at Re=40

圖6 Re=40，升力功率譜密度Fig.6 Power spectral density of lift at Re=40

在Re=200的工況下，Gera等[22]計(jì)算得到的阻力系數(shù)為1.487，文中模擬計(jì)算得到的阻力系數(shù)平均值為1.496，誤差為0.605%。圖7為Re=200 時(shí)計(jì)算得到的阻力系數(shù)和升力系數(shù)時(shí)歷曲線。在流動(dòng)充分發(fā)展之后，圖中升力系數(shù)曲線有明顯的周期性，升力變化頻率即為泄渦頻率。圖8為升力功率譜密度，最高的尖峰對(duì)應(yīng)的橫坐標(biāo)為斯特勞哈爾數(shù)St=0.14，與Okajima[3]試驗(yàn)結(jié)果0.145相比，誤差為3.57%。

圖7 Re=200，阻力系數(shù)和升力系數(shù)時(shí)歷曲線Fig.7 Time history of the drag and lift at Re=200

圖8 Re=200，升力功率譜密度Fig.8 Power spectral density of lift at Re=200

Re=1 000時(shí)，阻力系數(shù)和升力系數(shù)時(shí)歷曲線如圖9、圖10所示，計(jì)算得到的阻力系數(shù)為2.396，比王建春[23]用FVM方法計(jì)算的結(jié)果大7.44%，斯特勞哈爾數(shù)St為0.124，與王建春[23]的計(jì)算結(jié)果及前人試驗(yàn)結(jié)果相近。

圖9 Re=1 000，阻力系數(shù)和升力系數(shù)時(shí)歷曲線Fig.9 Time history of the drag and lift at Re=1 000

圖10 Re=1 000，升力功率譜密度Fig.10 Power spectral density of lift at Re=1 000

3.2 尾流流態(tài)分析

Re=40時(shí)，流場(chǎng)充分發(fā)展后，流動(dòng)會(huì)趨于穩(wěn)定，方柱后面出現(xiàn)一對(duì)上下對(duì)稱的尾渦，不會(huì)發(fā)生渦脫落等現(xiàn)象，與前人試驗(yàn)及計(jì)算結(jié)果一致，圖11為流場(chǎng)穩(wěn)定后方柱周圍流場(chǎng)示意。從圖11可以看到，方柱周圍流場(chǎng)是對(duì)稱的，方柱前后速度較低，兩側(cè)速度較大。圖12為模擬計(jì)算得到流場(chǎng)穩(wěn)定后流線圖，可以明顯觀察到方柱后方由于回流形成的對(duì)稱尾渦，與Sen等[21]在文獻(xiàn)中給出的結(jié)果相似。圖13為y/D=6處即方柱中心線上x方向時(shí)均速度，圖中速度小于0區(qū)域即為回流區(qū)，Re=40時(shí)回流區(qū)長(zhǎng)度L與方柱邊長(zhǎng)D的比值L/D=2.38，根據(jù)Sen等[21]給出的經(jīng)驗(yàn)公式L/D=-0.078 3+0.072 4×Re(5

圖11 Re=40方柱周圍流場(chǎng)Fig.11 Flow field around the square at Re=40

圖12 Re=40方柱周圍流線對(duì)比Fig.12 Stream-line around the square at Re=40

圖13 不同雷諾數(shù)下y/D=6處x方向時(shí)均速度Fig.13 Time-mean velocity in x direction of y/D=6

圖14 Re=40時(shí)，不同流向位置處x方向時(shí)均速度Fig.14 Velocity in x direction at different positions at Re=40

Re=200工況下，流場(chǎng)充分發(fā)展后，方柱后會(huì)出現(xiàn)旋渦周期性生成、脫落和卡門渦街現(xiàn)象。圖15為兩個(gè)不同時(shí)刻方柱周圍速度場(chǎng)。從圖15中可以看到，隨著來流速度的增加，方柱周圍流場(chǎng)越來越復(fù)雜，在雷諾數(shù)200時(shí)，方柱周圍速度場(chǎng)不再呈現(xiàn)對(duì)稱的性質(zhì)，而是尾跡不停擺動(dòng)。

圖15 Re=200，方柱周圍流場(chǎng)Fig.15 The flow field around the square at Re=200

圖16展示了Re=200時(shí)一個(gè)周期里不同時(shí)刻方柱周圍的流線，圖17為不同時(shí)刻的渦量云圖。從圖16可以明顯看到一個(gè)周期內(nèi)，方柱尾渦的生成和脫落，一個(gè)周期開始時(shí)，方柱左下角的渦即將脫落；T/5時(shí)右上角的渦不斷發(fā)展壯大，同時(shí)左下角重新生成小渦；之后左下角的小渦發(fā)展壯大，右上角的大渦開始遷移脫落，到2T/5時(shí)，右上角的渦完全脫落；然后左下角的渦遷移脫落，右上角生成小渦，再發(fā)展壯大，T時(shí)刻的流線與開始時(shí)刻幾乎相同，方柱尾流中的旋渦完成一個(gè)周期的運(yùn)動(dòng)。此后尾渦將重復(fù)上面過程，不斷生成小渦，發(fā)展，遷移，脫落，在下游形成著名的卡門渦街現(xiàn)象，從圖17渦量云圖中可以清晰地看到卡門渦街現(xiàn)象。

圖16 Re=200，一個(gè)周期內(nèi)方柱周圍流線Fig.16 Stream-line around the square in a period at Re=200

圖17 Re=200，整體渦量云圖Fig.17 Voticity contours at Re=200

在MPS方法中，粒子按拉格朗日描述法進(jìn)行自由運(yùn)動(dòng)，可以跟蹤粒子每個(gè)時(shí)間步的位置，得到粒子的運(yùn)動(dòng)軌跡，即跡線圖。圖18為Re=200時(shí)，某兩個(gè)粒子的跡線，從跡線中可以明顯看到由于繞流的發(fā)生，當(dāng)粒子運(yùn)動(dòng)到方柱后方時(shí)，發(fā)生了回流的現(xiàn)象。

當(dāng)Re=1 000時(shí)，與Re=200時(shí)觀察到的方柱繞流現(xiàn)象相似，方柱后方仍存在渦交替脫落現(xiàn)象。圖19為方柱周圍平均壓強(qiáng)系數(shù)CP與文獻(xiàn)[24]結(jié)果對(duì)比，A-B與B-C兩個(gè)區(qū)域的平均壓強(qiáng)系數(shù)與前人模擬結(jié)果吻合較好，C-D區(qū)域結(jié)果比文獻(xiàn)結(jié)果略大。圖20為方柱周圍渦量圖與王建春[23]的模擬結(jié)果對(duì)比，對(duì)比結(jié)果較好。還可以看到，由于流速增加，對(duì)流影響加強(qiáng)，渦的穩(wěn)定性受到影響，兩個(gè)渦脫周期內(nèi)的渦脫落形式有一定的差別。圖13中可以看到雷諾數(shù)為1 000時(shí)的回流區(qū)長(zhǎng)度小于雷諾數(shù)為200時(shí)的長(zhǎng)度，可能也是因?yàn)闇u穩(wěn)定性降低，脫落需要的發(fā)展長(zhǎng)度變短。

圖19 Re=1 000，方柱周圍平均壓強(qiáng)系數(shù)Fig.19 Mean-pressure coeffient over the cylinder at Re=1 000

圖20 Re=1 000時(shí)方柱周圍渦量云圖Fig.20 Voticity contours around the square at Re=1 000

4 結(jié) 語

運(yùn)用MPS方法，建立入口、出口邊界條件，精確模擬出均勻流。在此基礎(chǔ)上，成功模擬了雷諾數(shù)分別為40、200、1 000的二維方柱繞流。計(jì)算結(jié)果表明：

1) 在雷諾數(shù)為40時(shí)，方柱后方的尾渦不會(huì)脫落，隨著雷諾數(shù)增加，方柱繞流流場(chǎng)流動(dòng)狀態(tài)由定常渦轉(zhuǎn)變?yōu)橹芷谛孕郎u脫落狀態(tài)，雷諾數(shù)為200和1 000時(shí)，尾渦都會(huì)脫落，出現(xiàn)明顯的卡門渦街，且雷諾數(shù)為1 000 時(shí)渦的穩(wěn)定性比雷諾數(shù)為200時(shí)低，回流區(qū)長(zhǎng)度比雷諾數(shù)為200時(shí)短。

2) 用MPS方法計(jì)算得到的方柱升阻力系數(shù)、斯特勞哈爾數(shù)等特性與前人試驗(yàn)結(jié)果及計(jì)算結(jié)果接近，流場(chǎng)特征與前人結(jié)果相符，驗(yàn)證了MPS方法模擬低雷諾數(shù)下方柱繞流問題的可靠性及有效性。

但從圖17和圖20渦量云圖中可以看到，模擬結(jié)果存在一定的數(shù)值噪聲，造成這種噪聲的原因之一可能是粒子分布不完全均勻。而且如果雷諾數(shù)繼續(xù)增大，流速增加會(huì)加劇這種分布不均，進(jìn)而造成模擬結(jié)果偏差較大。引入粒子位移修正技術(shù)(particle shifting technique)可能可以解決這個(gè)問題，減小模擬誤差的同時(shí)增大可模擬的雷諾數(shù)。