基于MPS方法模擬低雷諾數方柱繞流問題

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

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

柱體繞流問題是流體力學領域的經典問題和研究熱點之一。當流體流過柱體時，會發生分離流、渦流等很多復雜的流動現象，準確地模擬這些現象一直是計算流體力學的目標之一。同時柱體繞流也是現代工程中經常遇到的問題，風吹過高層建筑物，水流流過橋墩，海流流過海洋平臺等都是典型的柱體繞流。當流體流過這些結構物時，結構物會受到很大的作用力，對結構的強度提出了更高的要求。

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

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

1 數值方法

1.1 控制方程

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

(1)

(2)

式中：P、V分別為壓力和速度，ν、ρ分別表示流體的運動黏性系數和密度，f表示質量力。MPS方法的控制方程為拉格朗日形式，因此不存在對流項。

1.2 核函數

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

(3)

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

1.3 不可壓縮條件

采用被Lee等[20]改寫過的混合源項法，表達式為：

(4)

(5)

其中，r是粒子的坐標，i、j為粒子的編號。

1.4 梯度模型

采用的梯度模型表達式為：

(6)

1.5 散度模型

在MPS方法的連續性方程中存在散度項，需要用核函數對其進行離散。采用的散度模型表達式為：

(7)

1.6 Laplacian模型

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

(8)

其中，λ是修正系數，引入λ是為了修正數值結果，使其與擴散方程的解析結果相一致，其表達式為：

(9)

2 數值模擬

2.1 計算模型

研究中，計算域示意如圖1所示。設定方柱邊長D為0.1 m，計算域大小設置為25D×12D。為消除入口邊界對柱體繞流模擬的影響，計算域流體入口距柱體中心7D，計算域出口與柱體中心的距離為18D，從而確保尾流得到充分發展。為減少壁面邊界對柱體周圍流場的影響，兩側壁面距柱體中心的距離均設置為6D。計算坐標系原點在計算域左下角，x軸正方向為來流速度方向。粒子間距設置為Δx/D=25，計算雷諾數分別為Re=40、200及1 000。

圖1 計算域示意Fig.1 Computational domain of flow around square

2.2 邊界條件

2.2.1 入口邊界條件

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

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

2.2.2 出口邊界條件

如圖2所示，當流體粒子流出計算域時，把流體粒子轉換成虛粒子，不影響其它流體粒子。虛粒子的質量、密度等物理量都是0，等待替補到入口邊界。虛粒子填入由于推板后退產生的空隙中后，會被賦予相應的物理性質，變成新的流體粒子，重新進入計算域。

2.2.3 壁面邊界條件

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

2.3 來流條件

文中來流條件設定為均勻來流，通過模擬管道內均勻流對均勻來流條件進行驗證。選取柱體繞流雷諾數Re=40的工況，設定來流速度為0.000 4 m/s，模擬沒有柱體時的管道流動。圖3為模擬得到的整體速度云圖，管道整體速度都在0.000 4 m/s附近。圖4為流動充分發展后x/D=7、9、12、17、22處的x方向時均速度，其中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 結果分析

3.1 水動力系數

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

(10)

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

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

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

圖5 Re=40，阻力系數和升力系數時歷曲線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]計算得到的阻力系數為1.487，文中模擬計算得到的阻力系數平均值為1.496，誤差為0.605%。圖7為Re=200 時計算得到的阻力系數和升力系數時歷曲線。在流動充分發展之后，圖中升力系數曲線有明顯的周期性，升力變化頻率即為泄渦頻率。圖8為升力功率譜密度，最高的尖峰對應的橫坐標為斯特勞哈爾數St=0.14，與Okajima[3]試驗結果0.145相比，誤差為3.57%。

圖7 Re=200，阻力系數和升力系數時歷曲線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時，阻力系數和升力系數時歷曲線如圖9、圖10所示，計算得到的阻力系數為2.396，比王建春[23]用FVM方法計算的結果大7.44%，斯特勞哈爾數St為0.124，與王建春[23]的計算結果及前人試驗結果相近。

圖9 Re=1 000，阻力系數和升力系數時歷曲線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 尾流流態分析

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

圖11 Re=40方柱周圍流場Fig.11 Flow field around the square at Re=40

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

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

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

Re=200工況下，流場充分發展后，方柱后會出現旋渦周期性生成、脫落和卡門渦街現象。圖15為兩個不同時刻方柱周圍速度場。從圖15中可以看到，隨著來流速度的增加，方柱周圍流場越來越復雜，在雷諾數200時，方柱周圍速度場不再呈現對稱的性質，而是尾跡不停擺動。

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

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

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

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

在MPS方法中，粒子按拉格朗日描述法進行自由運動，可以跟蹤粒子每個時間步的位置，得到粒子的運動軌跡，即跡線圖。圖18為Re=200時，某兩個粒子的跡線，從跡線中可以明顯看到由于繞流的發生，當粒子運動到方柱后方時，發生了回流的現象。

當Re=1 000時，與Re=200時觀察到的方柱繞流現象相似，方柱后方仍存在渦交替脫落現象。圖19為方柱周圍平均壓強系數CP與文獻[24]結果對比，A-B與B-C兩個區域的平均壓強系數與前人模擬結果吻合較好，C-D區域結果比文獻結果略大。圖20為方柱周圍渦量圖與王建春[23]的模擬結果對比，對比結果較好。還可以看到，由于流速增加，對流影響加強，渦的穩定性受到影響，兩個渦脫周期內的渦脫落形式有一定的差別。圖13中可以看到雷諾數為1 000時的回流區長度小于雷諾數為200時的長度，可能也是因為渦穩定性降低，脫落需要的發展長度變短。

圖19 Re=1 000，方柱周圍平均壓強系數Fig.19 Mean-pressure coeffient over the cylinder at Re=1 000

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

4 結 語

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

1) 在雷諾數為40時，方柱后方的尾渦不會脫落，隨著雷諾數增加，方柱繞流流場流動狀態由定常渦轉變為周期性旋渦脫落狀態，雷諾數為200和1 000時，尾渦都會脫落，出現明顯的卡門渦街，且雷諾數為1 000 時渦的穩定性比雷諾數為200時低，回流區長度比雷諾數為200時短。

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

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