用移動粒子半隱式方法數值模擬Poiseuille 流動問題

蘭小杰，趙偉文，萬德成

1 上海交通大學 船海計算水動力學研究中心，上海 200240

2 上海交通大學 船舶海洋與建筑工程學院，上海 200240

0 引 言

Poiseuille 流動問題廣泛存在于工業生產中，其最典型的例子是低雷諾數時的管道流動。近年來，隨著在深海資源開采方面的競爭越來越激烈，礦物管道輸運系統作為深海采礦的關鍵環節，對其的研究需求愈發迫切。而針對Poiseuille 流動的機理研究對于更復雜的管道運輸問題來說具有一定的指導意義。

早在19 世紀，Hagen 和Poiseuille 就從實驗中歸納出了低雷諾數圓管中的液體層性管流規律，即圓管截面上的速度分布為拋物線分布[1]。1968 年，Fox 等[2]通過實驗對Poiseuille 流的穩定性進行了研究，發現當雷諾數Re＞ 2 150 時，Poiseuille 流開始出現失穩現象。Papanastasiou 等[3]使用分離變量的方法求得了Poiseuille 流問題的理論解。陳雷等[4]對不同邊界條件下的非穩態不可壓縮Poiseuille 流的發展過程進行了研究，發現對于平均速度隨著時間從0 線性增加邊界、恒壓力邊界、恒平均速度邊界這3 種不同類型的邊界條件，對應的非穩態發展過程依次縮短。金開文等[5]使用格子Boltzmann 方法對Poiseuille 流進行了模擬研究，發現模擬結果與理論結果吻合，說明采用格子Boltzmann 方法處理壓力驅動類層流問題具有可行性。在粒子法方面，Adami 等[6]采用光滑粒子流體動力學方法（ SPH）對Re=0.001 25 時的二維Poiseuille流進行了模擬，發現模擬結果與理論結果吻合較好，證明采用SPH 方法模擬低雷諾數下的Poiseuille 流是有效的。Meister 等[7]利用弱可壓縮SPH 方法對Poiseuille 流進行了模擬，發現將SPH 方法應用于中、高雷諾數下的Poiseuille 流（Re≥ 1）會導致橫向不穩定性問題，并對造成這種不穩定性的原因進行了探討。劉謀斌和常建忠[8]采用SPH 方法模擬了Poiseuille 流，發現其穩定一段時間后會出現模擬結果逐漸偏離穩定解的情況，說明采用原有SPH 方法模擬Poiseuille 流會存在數值不穩定的情況，然后采用一種改進的SPH 方法— 有限粒子法模擬了Poiseuille 流，結果發現數值不穩定的情況消失了，可見，有限粒子法是解決原有SPH 方法數值不穩定的一個有效方法。Song 等[9]利用SPH 方法對Re= 0.01～100 下的Poiseuille 流進行了模擬，系統地研究了參數、背景壓力、初始粒子密度和密度重新初始化技術對SPH 模擬Poiseuille 流的影響。綜上所述，粒子法對Poiseuille 流的研究大多集中在SPH 方法上，鮮有人采用移動粒子半隱式（moving particle semi-implicit，MPS）方法對Poiseuille流進行研究。

1996 年，Koshizuka 和Oka[10]在SPH 方法的基礎上提出了MPS 方法，主要用于模擬不可壓縮流動問題。與網格類方法最大的不同之處是，在MPS 方法中，流體是通過相互作用的粒子來離散的，并利用拉格朗日粒子攜帶空間流場的信息，粒子之間的影響則通過核函數來實現。由于粒子運動不會受到網格間固定拓撲關系的限制，因此在處理大變形的自由面問題時不會出現網格畸變的問題，具有更大的靈活性。粒子法目前已被廣泛應用于自由表面的流動，以及多相流和水下爆炸等流動中[11-16]，但將其應用于壁面剪切流動時的數值穩定性還需進一步予以驗證。為此，本文將使用MPS 求解器MLParticle-SJTU，通過建立恒流量入口邊界條件和無滑移壁面邊界條件，模擬二維Poiseuille 流并與理論解析解進行對比，驗證MPS 方法在模擬Poiseuille 流問題上的有效性。

1 數值方法

1.1 控制方程

在MPS 方法中，控制方程一般包括連續性方程以及N-S 方程，分別可以寫成如下形式：

式中：ρ 為流體密度，m3/kg；V為流體速度矢量，m/s；P為流體壓力，Pa； ν為流體運動黏性系數，m2/s；f為流體質量力，m/s2。

1.2 核函數

有別于SPH 方法，MPS 方法中的核函數只是作為權函數來使用，其求解過程無需使用核函數的導函數，所以MPS 方法的核函數只要求連續而不要求光滑。本研究采用的核函數為文獻[15]推薦的核函數，表達式如下：

1.3 不可壓縮條件

本文采用Lee 等[17]改寫的混合源項法，表達式如下：

式中：γ 為泊松方程源項中粒子數密度的權重，可取0～1 之間的任意數值，本文取γ = 0.99；為粒子i的臨時速度矢量； 在MPS 方法中，通常使用粒子數密度來表示粒子的疏密，n0為初始粒子數密度，n*為臨時粒子數密度。＜n*＞的表達式如下：

1.4 梯度模型

本文采用的梯度模型為：

式中： φ為任意標量， φi和 φj為此標量在粒子i，j處的值；d為計算維數，本文研究的是二維Poiseuille 流。

1.5 散度模型

在連續性方程中存在散度項，需使用核函數對散度項進行離散。本文采用的散度模型為：

式中， Π為任意矢量， Πi和 Πj分別為粒子i，j處該矢量的值。

1.6 Laplacian 模型

本文使用Koshizuka 和Oka[10]所給的如下Laplacian 模型：

式中，λ 為修正系數。引入λ 可以修正數值計算的結果，使其與擴散方程的解析結果一致，其表達式為：

1.7 邊界條件

如圖1 所示，入口采取推板的形式形成恒流量入口，推板向前推動流體流動，在到達指定的位置后，推板退回至初始位置，再在空出的位置填入幽靈粒子，并賦予它們質量、黏性、密度等與流體相關的物理量，使其轉換成流體粒子。將從出口流出計算域的流體粒子轉換成幽靈粒子，在入口推板后退時填入入口處，并重新轉換回流體粒子，幽靈粒子的質量、密度等物理量均為0。MPS 方法中的壁面邊界由多層粒子組成，與流體顆粒相鄰的邊界粒子為第1 類邊界粒子，其壓力通過壓力Poisson 方程得到；不與流體粒子接觸的粒子為第2 類邊界粒子，其壓力通過周圍流體粒子和第1 類邊界粒子向外插值得到。

圖1 邊界示意圖Fig. 1 Schematic diagram of boundaries

上、下邊界均設置為不可滑移壁面，圖2 所示為無滑移壁面示意圖。圖中：Ui，u⊥，u||分別為流體粒子的速度、垂直于壁面的速度和平行于壁面的速度；′分別為與流體粒子對應的虛擬粒子的速度、垂直于壁面的速度和平行于壁面的速度，這些速度的單位均為m/s。在每個時間步內計算黏性力時，在壁面外側生成虛擬粒子，虛擬粒子與流體粒子關于壁面軸對稱。在無滑移壁面條件下，當壁面靜止時，虛擬粒子的垂向速度和切向速度均與流體粒子相反；在計算黏性力時，需考慮這些虛擬粒子對流體粒子的作用力，從而形成無滑移壁面邊界條件。

圖2 不可滑移邊界示意圖Fig. 2 Schematic diagram of no-slip boundary

2 數值模擬

2.1 Poiseuille 流

所謂Poiseuille 流，是指由壓力驅動的層性管流，在流動充分發展后，截面上的速度分布將呈拋物線形狀，如圖3 所示。圖中，D為管道直徑。

圖3 Poiseuille 流示意圖Fig. 3 Schematic diagram of Poiseuille flow

由于Poiseuille 流動條件簡單，其基本方程組可以有解析解，解析解如下：

式中：u為流體在x方向的速度；Q為流量。由式（10）可以看出，Poiseuille 流的速度分布呈拋物線形狀，通過流速的分布，可以驗證后續經MPS 方法模擬得到的結果。

2.2 計算模型

圖4 所示為計算域的示意圖。設D= 0.2 m，為了讓流動能夠充分發展，將計算流體域的長度設置為15D。此外，設流體的密度ρ= 1 000 kg/m3，黏性系數μ= 0.001 Pa·s。計算工況包括3 個，即入口速度U= 0.000 02，0.000 2 和0.002 m/s，其雷諾數Re= 4，40，400。

圖4 計算域示意圖Fig. 4 Schematic diagram of computational domain

2.3 收斂性分析

選取Re= 400 的工況，采用3 種不同的粒子間距dx= 0.005 5，0.004 和0.002 5 m 來對粒子的收斂性進行驗證，3 種粒子間距對應的流體粒子數量分別為18 410，36 750 和94 800。圖5 對不同粒子間距下流動充分發展后的截面速度分布與解析解進行了對比，圖中，y/D為截面縱坐標與管道直徑的比值。由圖可見，不同粒子間距對應的截面最大速度分別為Vx= 0.002 62，0.002 875，0.002 822 m/s，與解析解的差距分別為12.7%，4.0%和5.9%；dx= 0.004 m 與dx= 0.002 5 m 的計算結果十分接近，說明當dx達到0.004 m 后，繼續減小粒子間距對模擬結果的影響很小，考慮到計算量，接下來的工況將均選取dx= 0.004 m 來進行相關計算。

圖5 Re = 400 時不同粒子間距下截面x 方向的速度分布Fig. 5 Velocity distribution of flow in the x direction in different dx when Re = 400

2.4 模擬結果

圖6 所示為不同雷諾數下，流動充分發展后流體域在x方向的速度云圖。從圖中可以看到，在無滑移壁面及流體黏性的作用下，3 個工況下的速度均呈現出上下界面的速度小、中間的速度最大的規律。圖7 所示為不同雷諾數下流動充分發展后中心線上的速度分布，圖中Vx/U為x方向速度與入口速度的比值。從中可以看到，隨著Re的增加，流動充分發展所需長度也隨之增加，即Re= 4，40，400 時流動充分發展所需長度分別約為1.5D，2.5D和10D。

圖6 流體域在x 方向的速度云圖Fig. 6 Velocity contours of fluid domain in the x direction

圖7 流動充分發展后在中心線上x 方向的速度分布Fig. 7 Velocity distribution of fully-developed flow on the central axis in the x direction

選取x=10D的截面，觀察截面處x方向的速度分布隨時間的變化，如圖8 所示。由圖可見，在不同工況下及流動發展的初始階段，截面不同位置處的流速差均較??；隨著流動的發展，中間位置的流速不斷增大；待流動充分發展后，速度剖面曲線最終呈拋物線形狀。在Re= 4 工況下，當無量綱時間Ut/D= 0.5 時，速度剖面曲線與理論解析解幾乎重合，此時流動發展到穩定狀態，至Ut/D= 10 時，截面的速度峰值有所回落，其值略小于理論解析解，相對誤差為2.6%，不過也只是在理論解析解附近有小幅度的振蕩，未有完全偏離穩定解的情況，說明采用MPS 方法模擬Poiseuille流動問題不存在數值不穩定現象。在Re= 40 工況下，當Ut/D= 2 時，截面處x方向的速度分布與Ut/D= 10 時的幾乎完全一致，說明流動已充分發展，此時截面處的最大速度Vmax= 0.000 293 4 m/s，略小于0.000 3 m/s 的理論解析解，相對誤差為2.1%。在Re= 400 時，Vmax= 0.002 875 m/s，與解析結果間的相對誤差為4.0%。

圖8 x = 10D 截面處x 方向的速度分布隨時間的變化Fig. 8 Velocity distribution of flow in the x direction changes with time at the section of x = 10D

在MPS 方法中，粒子按拉格朗日描述法進行自由運動，其可以跟蹤粒子每個時間步的位置。選取Re= 4 的工況，跟蹤初始時刻x方向位置在x= 10D截面處流體粒子的運動，圖9 所示為其在不同時刻的分布形狀。從中可以看到，在初始時刻，同一截面處的粒子由于運動速度不一致，隨著時間的發展會呈現拋物線形狀。

圖9 Re = 4 工況下在初始時刻x = 10D 截面處粒子隨時間的變化分布Fig. 9 Distribution of particle moves with time at the section of x =10D at initial time when Re = 4

3 結 語

本文運用MPS 方法，建立了出入口邊界及無滑移壁面條件，并對不同雷諾數下的二維Poiseuille流進行了模擬。計算結果表明：Poiseuille 流充分發展后，流體速度的分布情況是從中間向兩側逐漸遞減，速度分布呈拋物線形狀，且在恒流量入口邊界條件下，隨著雷諾數的升高，Poiseuille 流發展所需要的流動長度變大；在不同雷諾數下，Poiseuille 流經充分發展后，模擬得到的速度最大值與理論解析解之間的相對誤差在5%以內，驗證了采用MPS 方法模擬Poiseuille 流問題的可靠性及有效性。下一步，將對采用MPS 方法模擬相對更復雜的柱體繞流相關問題的可靠性進行計算與驗證。