用移動(dòng)粒子半隱式方法數(shù)值模擬Poiseuille 流動(dòng)問(wèn)題

蘭小杰，趙偉文，萬(wàn)德成

1 上海交通大學(xué) 船海計(jì)算水動(dòng)力學(xué)研究中心，上海 200240

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

0 引 言

Poiseuille 流動(dòng)問(wèn)題廣泛存在于工業(yè)生產(chǎn)中，其最典型的例子是低雷諾數(shù)時(shí)的管道流動(dòng)。近年來(lái)，隨著在深海資源開(kāi)采方面的競(jìng)爭(zhēng)越來(lái)越激烈，礦物管道輸運(yùn)系統(tǒng)作為深海采礦的關(guān)鍵環(huán)節(jié)，對(duì)其的研究需求愈發(fā)迫切。而針對(duì)Poiseuille 流動(dòng)的機(jī)理研究對(duì)于更復(fù)雜的管道運(yùn)輸問(wèn)題來(lái)說(shuō)具有一定的指導(dǎo)意義。

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

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

1 數(shù)值方法

1.1 控制方程

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

式中：ρ 為流體密度，m3/kg；V為流體速度矢量，m/s；P為流體壓力，Pa； ν為流體運(yùn)動(dòng)黏性系數(shù)，m2/s；f為流體質(zhì)量力，m/s2。

1.2 核函數(shù)

有別于SPH 方法，MPS 方法中的核函數(shù)只是作為權(quán)函數(shù)來(lái)使用，其求解過(guò)程無(wú)需使用核函數(shù)的導(dǎo)函數(shù)，所以MPS 方法的核函數(shù)只要求連續(xù)而不要求光滑。本研究采用的核函數(shù)為文獻(xiàn)[15]推薦的核函數(shù)，表達(dá)式如下：

1.3 不可壓縮條件

本文采用Lee 等[17]改寫(xiě)的混合源項(xiàng)法，表達(dá)式如下：

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

1.4 梯度模型

本文采用的梯度模型為：

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

1.5 散度模型

在連續(xù)性方程中存在散度項(xiàng)，需使用核函數(shù)對(duì)散度項(xiàng)進(jìn)行離散。本文采用的散度模型為：

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

1.6 Laplacian 模型

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

式中，λ 為修正系數(shù)。引入λ 可以修正數(shù)值計(jì)算的結(jié)果，使其與擴(kuò)散方程的解析結(jié)果一致，其表達(dá)式為：

1.7 邊界條件

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

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

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

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

2 數(shù)值模擬

2.1 Poiseuille 流

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

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

由于Poiseuille 流動(dòng)條件簡(jiǎn)單，其基本方程組可以有解析解，解析解如下：

式中：u為流體在x方向的速度；Q為流量。由式（10）可以看出，Poiseuille 流的速度分布呈拋物線形狀，通過(guò)流速的分布，可以驗(yàn)證后續(xù)經(jīng)MPS 方法模擬得到的結(jié)果。

2.2 計(jì)算模型

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

圖4 計(jì)算域示意圖Fig. 4 Schematic diagram of computational domain

2.3 收斂性分析

選取Re= 400 的工況，采用3 種不同的粒子間距dx= 0.005 5，0.004 和0.002 5 m 來(lái)對(duì)粒子的收斂性進(jìn)行驗(yàn)證，3 種粒子間距對(duì)應(yīng)的流體粒子數(shù)量分別為18 410，36 750 和94 800。圖5 對(duì)不同粒子間距下流動(dòng)充分發(fā)展后的截面速度分布與解析解進(jìn)行了對(duì)比，圖中，y/D為截面縱坐標(biāo)與管道直徑的比值。由圖可見(jiàn)，不同粒子間距對(duì)應(yīng)的截面最大速度分別為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 的計(jì)算結(jié)果十分接近，說(shuō)明當(dāng)dx達(dá)到0.004 m 后，繼續(xù)減小粒子間距對(duì)模擬結(jié)果的影響很小，考慮到計(jì)算量，接下來(lái)的工況將均選取dx= 0.004 m 來(lái)進(jìn)行相關(guān)計(jì)算。

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

2.4 模擬結(jié)果

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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