多相流動(dòng)的光滑粒子流體動(dòng)力學(xué)方法研究綜述1)

陳飛國(guó) * 葛 蔚

* (中國(guó)科學(xué)院過(guò)程工程研究所多相復(fù)雜系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100190)

? (中國(guó)科學(xué)院綠色過(guò)程制造創(chuàng)新研究院,北京 100190)

** (中國(guó)科學(xué)院大學(xué)化學(xué)工程學(xué)院,北京 100049)

引言

多相流廣泛存在于自然界和工業(yè)應(yīng)用中.對(duì)多相流進(jìn)行高精度和高效率的建模和模擬非常重要.然而,其中流動(dòng)形態(tài)和界面跟蹤的復(fù)雜性限制了基于網(wǎng)格的傳統(tǒng)方法模擬多相流能力,尤其對(duì)于存在界面大變形、不連續(xù)性和多物理過(guò)程的多相流動(dòng).而無(wú)網(wǎng)格方法可以很好地處理網(wǎng)格方法遇到的上述問(wèn)題.在過(guò)去幾十年發(fā)展起來(lái)的無(wú)網(wǎng)格方法中,光滑粒子流體動(dòng)力學(xué)方法(smoothed particle hydrodynamics,SPH)是眾所周知并被廣泛應(yīng)用的其中一種.

SPH 方法由Gingold 和Monaghan[1]以及Lucy[2]于1977 年發(fā)明,最初用于模擬天體物理學(xué)中的非軸對(duì)稱(chēng)現(xiàn)象.它已獲得極大發(fā)展并成為模擬流體運(yùn)動(dòng)的標(biāo)準(zhǔn)方法之一[3].已有包含各種綜述的大量文獻(xiàn)[3-20](相關(guān)文獻(xiàn)的數(shù)量極其巨大,Google Scholar 可查詢(xún)到近10 萬(wàn)個(gè)結(jié)果)論述了SPH 方法各方面及在各領(lǐng)域中的應(yīng)用情況,如界面流動(dòng)、流變學(xué)、流動(dòng)不穩(wěn)定性、多相流、破損、彈性和斷裂力學(xué)、爆炸力學(xué)、波傳播和熱傳導(dǎo)等.

從純粹的拉格朗日法角度來(lái)看,SPH 方法具有以下主要優(yōu)點(diǎn)[7,14]:

1)易于處理大變形和快速移動(dòng)的自由表面;

2)自發(fā)捕獲和跟蹤間斷面,而這是使用基于網(wǎng)格的方法難以實(shí)現(xiàn)的;

3)能有效模擬耦合多物理場(chǎng)的多相流;

4) SPH 粒子大小可靈活調(diào)節(jié),對(duì)應(yīng)不同的數(shù)值精度要求;

5)二維和三維SPH 模擬均易于實(shí)現(xiàn)并具有良好的并行計(jì)算可擴(kuò)展性.

上述特性確保了SPH 方法特別適用于基于網(wǎng)格的方法很難或者無(wú)法完成的多相流高精度模擬.許多研究人員已應(yīng)用 SPH 方法模擬多相流,這些應(yīng)用中SPH 的實(shí)現(xiàn)方式有很大不同.最常見(jiàn)的是多相SPH 方法[21-29],其中兩種流體分別采用SPH 方法求解.此外,SPH 方法還可被視為拉格朗日求解器求解表征兩相流動(dòng)的雙流體模型(two-fluid model,TFM)[30-32].SPH 方法也可以與其他離散方法結(jié)合使用,多用于含固體顆粒的多相流動(dòng)模擬.為了降低計(jì)算成本,基于網(wǎng)格的歐拉求解器耦合到 SPH 方法中,處理不存在相界面的流動(dòng)主體部分.一些研究人員還發(fā)展了各種修正方法來(lái)提高SPH 方法在多相流模擬中的準(zhǔn)確性和穩(wěn)定性.

盡管眾多文獻(xiàn)涵蓋了 SPH 方法模擬各種多相流動(dòng),也有一些綜述文章論述了SPH 方法用于多相流動(dòng)模擬的各類(lèi)模型,但大多僅涉及多相流模擬的某類(lèi) SPH 方法實(shí)現(xiàn),很少總結(jié)SPH 方法在多相流模擬中的不同實(shí)現(xiàn)方式.本文將歸納和比較這些不同SPH 方法實(shí)現(xiàn)方式,并對(duì)相應(yīng)適用體系進(jìn)行建議.

本文組織如下:第1 節(jié)簡(jiǎn)單介紹SPH 方法;第2 節(jié)介紹模擬多相流的其他SPH 類(lèi)方法;第3 節(jié)總結(jié)和比較SPH 方法在多相流模擬中的各種實(shí)現(xiàn)方式;第4 節(jié)論述SPH 方法的各種改進(jìn)措施;最后在第5 節(jié)總結(jié)全文,并展望提高SPH 方法模擬多相流精度和效率的多種方式.

1 SPH 方法

SPH 方法是一種典型的拉格朗日粒子方法,它用帶有特定物理性質(zhì)(例如密度、黏度、質(zhì)量、速度和溫度)的粒子集合來(lái)表示流體[8,10,15].粒子的運(yùn)動(dòng)和性質(zhì)變化由質(zhì)量、動(dòng)量和能量守恒方程控制.它們的數(shù)值求解則由核近似和粒子插值的權(quán)重積分獲得.

1.1 SPH 方法的基本原理

簡(jiǎn)而言之,SPH 的基本思想[8]可概括如下:用一組帶有物理性質(zhì)的粒子來(lái)離散化問(wèn)題域,并將流體的控制方程轉(zhuǎn)變?yōu)榫哂懈鞣N邊界條件的常微分方程;然后,采用核函數(shù)插值近似的權(quán)重積分獲得數(shù)值解.SPH 方法的關(guān)鍵是插值近似,包括核近似和粒子近似.

1.1.1 核近似

問(wèn)題域Ω中點(diǎn)x的連續(xù)光滑函數(shù)值f(x)定義為[8]

為了消除由δ函數(shù)引入的無(wú)限值,引入光滑核函數(shù)來(lái)逼近δ函數(shù),則該函數(shù)值表達(dá)為[8,12]

其中W(x?x?,h)為核函數(shù),光滑長(zhǎng)度h用于衡量權(quán)重積分的影響區(qū)域.

函數(shù)導(dǎo)數(shù)由核函數(shù)近似表達(dá)為

1.1.2 粒子近似

問(wèn)題域Ω被離散化為一組粒子(如圖1 所示),粒子j處代表的空間體積即粒子j的體積為ΔVj.若粒子j的密度為ρj,則其質(zhì)量mj為

圖1 SPH 粒子近似Fig.1 Schematic illustration of SPH particle approximation

其在粒子i處的值為

這樣,函數(shù)的取值及其導(dǎo)數(shù)值可以通過(guò)離散粒子的加權(quán)平均得到.

1.2 控制方程的離散化

流體運(yùn)動(dòng)由守恒方程控制,即連續(xù)方程、動(dòng)量方程(納維?斯托克斯方程)和能量方程,它們以下列形式表示

其中u為速度,ρ為密度,S為應(yīng)力張量,g為單位質(zhì)量體積力,e為比內(nèi)能,q為熱通量.應(yīng)力張量S采用顯式形式

這里I表示單位張量,p為內(nèi)壓,剪切應(yīng)力σ表達(dá)為

其中μ和ξ分別為剪切黏度系數(shù)和體積黏度系數(shù).

守恒方程的SPH 表達(dá)式為[8-10,33]

在SPH 方法中,計(jì)算粒子密度有兩種形式,其中一種利用了密度的導(dǎo)數(shù)

需要注意的是,由于動(dòng)量方程中的壓力表達(dá)式未知,連續(xù)性方程和動(dòng)量方程并不封閉.流體的壓力一般由狀態(tài)方程(equation of state,EOS)獲得.對(duì)于不可壓縮SPH,則由投影方法求解壓力泊松方程來(lái)耦合速度和壓力.

1.3 狀態(tài)方程

在SPH 方法中,流體的狀態(tài)方程是有效計(jì)算動(dòng)量方程中壓力項(xiàng)的重要因素.早期的SPH 模擬[34]使用了 Batchelor 狀態(tài)方程[35]

其中γ是一個(gè)常數(shù),取值范圍一般從1 到7;下標(biāo)0 表示參考量.它也是弱可壓縮SPH 中最常用的狀態(tài)方程.對(duì)于較大的γ,壓力對(duì)密度的高敏感性可能導(dǎo)致密度波動(dòng)較大時(shí)壓力估算偏差極大.

另一個(gè)常用狀態(tài)方程則采用如下簡(jiǎn)單形式[36-37]

其中c是流體聲速;它是平衡真實(shí)流體近似值和時(shí)間演化效率的重要因素.

對(duì)于范德華(vdW)流體,狀態(tài)方程為[38]

其中n是數(shù)密度,T是溫度,kB是玻爾茲曼常數(shù),a和b分別是表達(dá)內(nèi)聚作用和分子大小的系數(shù);其中m為粒子質(zhì)量.狀態(tài)方程第二項(xiàng)是形成穩(wěn)定液滴表面的內(nèi)聚作用.范德華狀態(tài)方程也被稱(chēng)為內(nèi)聚壓力模型.

1.4 核函數(shù)

SPH 方法的核心問(wèn)題之一是如何有效地執(zhí)行準(zhǔn)確和穩(wěn)定的近似操作.用于核近似和粒子近似的光滑核函數(shù)決定了插值方式和加權(quán)積分.

許多研究人員對(duì)光滑核進(jìn)行了研究,并提出了多種核函數(shù).光滑核函數(shù)一般需要滿(mǎn)足如下要求[3,8-9,15]:

(1)在支持域上滿(mǎn)足正則化條件(歸一性)

(2)緊致性

式中rc表示光滑核函數(shù)有效作用距離.

(3)非負(fù)性

(4)當(dāng)光滑長(zhǎng)度趨于零時(shí),光滑核函數(shù)滿(mǎn)足狄拉克δ函數(shù)(δ函數(shù)性)

(5)遞減、偶對(duì)稱(chēng)以及充分光滑性質(zhì).

理論上,任何滿(mǎn)足上述要求的函數(shù)都可以作為SPH 方法的光滑核函數(shù).Liu 等[9,39]討論了它的一般構(gòu)造方法.一個(gè)常用的核函數(shù)是最初由 Gingold 和Monaghan[1]提出的高斯函數(shù)

使用最廣泛的是分段樣條函數(shù),如三次樣條函數(shù)[36,40]

其中s=|x?x?|/h,αd是與維度d相關(guān)的歸一化因子.有關(guān)典型光滑核函數(shù)的更多細(xì)節(jié)可參見(jiàn)已有的綜述文獻(xiàn)[8-9,15].

1.5 表面張力

存在于兩相界面處的表面張力是多相流的重要特性,特別是對(duì)于氣液流、微流動(dòng)流和表面流,它是影響流動(dòng)行為的主要因素.表面張力的準(zhǔn)確數(shù)值表達(dá)對(duì)于多相流動(dòng)建模至關(guān)重要.從分子相互作用的角度來(lái)看,表面張力源于界面上分子間力的不對(duì)稱(chēng)性[38],見(jiàn)圖2 所示.

圖2 分子在界面上的相互作用Fig.2 Schematic illustration of the interaction of molecules at the interface

在SPH 方法中,主要有三類(lèi)模型描述表面張力,即連續(xù)表面力(continuum surface force,CSF)模型、粒子間力(inter-particle force,IPF)模型和內(nèi)聚壓力模型(cohesive pressure model,CPM).有時(shí),CPM 被視為廣義的IPF 模型.

1.5.1 CSF 模型

CSF 模型由Brackbill 等[41]首先提出,它是一種宏觀(guān)力模型,表面張力被描述為兩相之間的偏向分子相互作用而兩相流體之間的分界面被視為有限厚度的過(guò)渡區(qū)域.表面張力的力密度定義為與表面的局部曲率成正比[37,42-45]并表達(dá)為

其中γ是表面張力系數(shù),κ 是界面的局部曲率,n是界面的單位法向量.

通常由下式定義的色函數(shù)區(qū)分不同的相[37,41]

然后將界面定義為具有非零色函數(shù)梯度的有限過(guò)渡帶.界面的法向量由下式計(jì)算

1.5.2 IPF 模型

IPF 模型中,表面張力由界面處SPH 粒子的對(duì)力隱式表達(dá).當(dāng)界面過(guò)渡區(qū)內(nèi)的兩個(gè)粒子i和j分別屬于不同相時(shí),IPF 通常表達(dá)為

與其他相粒子相互作用的附加力將導(dǎo)致由表面張力產(chǎn)生的凈加速度

文獻(xiàn)中常用的IPF 模型包括Lennard?Jones 力模型[46]

更多的IPF 模型見(jiàn)于Wang 等[15]的綜述.

1.5.3 CPM

CPM 可以被視為廣義 IPF 模型,其中界面粒子上的附加力由內(nèi)聚的范德華狀態(tài)方程隱式確定.假設(shè)內(nèi)聚壓力的相互作用范圍超過(guò)其他力的相互作用范圍[48];因此,內(nèi)部的粒子內(nèi)聚壓力項(xiàng)被平衡導(dǎo)致凈力為零,而界面區(qū)域的粒子內(nèi)聚壓力項(xiàng)不平衡由范德華狀態(tài)方程中的吸引項(xiàng)產(chǎn)生表面張力.界面區(qū)域內(nèi)粒子的凈加速度為[33,48]

1.6 固體邊壁

與其他離散方法一樣,SPH 方法中固壁邊界條件的實(shí)現(xiàn)并不像基于網(wǎng)格的方法那樣直接.SPH 方法中,實(shí)施固壁邊界條件的目的是防止流體粒子穿透固壁,并在光滑近似被截?cái)鄷r(shí)在邊界處進(jìn)行加權(quán)插值,從而實(shí)現(xiàn)自由滑移、非滑移或部分滑移邊界條件.離散方法(包括SPH 方法)中的固壁邊界處理方式原則上可以大致分為兩類(lèi),即粒子表示法和幾何表示法.

1.6.1 粒子表示法

在粒子表示法中,固體邊壁由固定粒子的集合表示,在流體粒子和邊界粒子之間施加附加力以產(chǎn)生指定的邊界條件.根據(jù)不同的粒子生成方式,固壁的粒子表示法又有邊界粒子模型[34]、鏡像粒子模型[49]和動(dòng)態(tài)粒子模型[8,50]等3 種典型模型.

在邊界粒子模型中,固體虛擬粒子位于邊界上,如圖3 所示.固體粒子施加到流體粒子上的排斥力用于防止流體非物理滲入固壁.研究人員提出了各種排斥力模型來(lái)改進(jìn)和修正邊界力計(jì)算[13,15,34].

圖3 邊界粒子模型Fig.3 Schematic illustration of solid virtual particles

鏡像粒子模型也被稱(chēng)為鬼粒子模型,其中鏡像粒子為流體粒子在固壁中的鏡像[49,51],如圖4 所示.鏡像粒子的位置和速度與流體粒子切向?qū)ΨQ(chēng),而其他屬性與流體粒子相同.因此,鏡像粒子信息會(huì)在每個(gè)時(shí)間步中動(dòng)態(tài)重新生成.

圖4 鏡像粒子模型Fig.4 Schematic illustration of image particles

在動(dòng)態(tài)粒子模型(圖5 所示)中,固體邊壁由邊界內(nèi)的一組固定粒子表示[8,50,52].邊界粒子的位置不隨時(shí)間變化而但其他屬性通過(guò)靠近邊界的流體粒子線(xiàn)性外推獲得.這些邊界內(nèi)粒子提供了動(dòng)態(tài)的邊界條件.

圖5 動(dòng)態(tài)粒子模型Fig.5 Schematic illustration of dynamics particles

1.6.2 幾何表示法

幾何表示法引入半解析模型[53-55].在半解析模型中(圖6 所示),邊界屬性(如線(xiàn)段和向內(nèi)法線(xiàn))由幾何代表性粒子重新定義[53].引入重整化因子來(lái)重建靠近邊界的流體粒子的控制方程.梯度算子則被表達(dá)為[53]

圖6 邊界屬性定義示意Fig.6 Schematic of boundary property definitions

2 多相流模擬其他SPH 類(lèi)方法

2.1 不可壓縮SPH

在常規(guī)SPH 方法中,由描述壓力和密度之間關(guān)系的狀態(tài)方程表達(dá)動(dòng)量方程中的壓力項(xiàng).因此,流體被認(rèn)為是可壓縮的而沒(méi)有速度無(wú)散約束.Cummins和Rudman[56]用投影方法求解壓力泊松方程,解決了SPH 方法中的速度無(wú)散問(wèn)題,從而建立不可壓縮SPH (incompressible SPH,ISPH)方法.

中間速度場(chǎng)u*是在沒(méi)有壓力梯度項(xiàng)的情況下在時(shí)間上向前積分動(dòng)量方程計(jì)算得到.壓力泊松方程基于中間速度場(chǎng)獲得并求解獲得達(dá)到不可壓縮條件的壓力場(chǎng)

Ellero 等[57]引入動(dòng)態(tài)約束使流體粒子體積恒定以實(shí)現(xiàn)不可壓縮性.Hu 和Adams[58]將投影方法成功應(yīng)用于多相流建立多相ISPH.Szewc 等[59]比較了ISPH 方法中基于網(wǎng)格的泊松求解器和基于粒子的泊松求解器.Cornelis 等[60]提出了隱式ISPH 用于水流的真實(shí)渲染.Aly[61]應(yīng)用ISPH 方法模擬了方管和方腔中的自然對(duì)流.Lind 等使用 ISPH 方法模擬了自由表面流動(dòng)[42]和空氣中液滴變形過(guò)程[62].

2.2 移動(dòng)粒子半隱式法

Koshizuka[63]針對(duì)核材料流體力學(xué)問(wèn)題提出了另一種基于拉格朗日的無(wú)網(wǎng)格方法——移動(dòng)粒子半隱式法(moving particle semi-implicit,MPS).隨著進(jìn)一步的發(fā)展,MPS 已被廣泛用于處理各類(lèi)問(wèn)題,如潰壩[64]、破波[65]、氣液兩相流[66-67]和自由表面流[68].Khayyer 和Gotoh[69]以及Kondo 和Koshizuka[70]修正MPS 中壓力項(xiàng),消除了非物理高頻數(shù)值振蕩.

MPS 具有與其他粒子方法相似的特征,控制方程中各項(xiàng)公式與SPH 方法非常相似.ISPH 和MPS都采用類(lèi)似的半隱式數(shù)值方法來(lái)求解不可壓縮流,它們之間的主要區(qū)別在于光滑核函數(shù)[71-72].在MPS中,經(jīng)常使用如下非常近程的光滑核函數(shù)[63]

它的值在r=0 時(shí)趨于無(wú)窮以避免粒子聚集,而ISPH中r=0 時(shí)的光滑核函數(shù)值為一有限值.

2.3 宏觀(guān)擬顆粒方法

宏觀(guān)擬顆粒方法(macro-scale pseudo-particle method,MaPPM) 是SPH 方法的另一種變體,由Ge 和Li[73-74]提出用于顆粒?流體系統(tǒng)高精度模擬.在MaPPM 中,各點(diǎn)處的導(dǎo)數(shù)計(jì)算被該點(diǎn)與其相鄰點(diǎn)的有限差分加權(quán)平均計(jì)算取代

式中,dim表示體系的維度,如是二維,dim=2;如是三維,則dim=3.

這種替換帶來(lái)了兩個(gè)主要好處:(1)在對(duì)多相流進(jìn)行模擬和分析時(shí),可以保持導(dǎo)數(shù)的簡(jiǎn)單加權(quán)平均性質(zhì);(2)可采用不可導(dǎo)函數(shù)作為光滑核函數(shù).

Zhou 等[29]在MaPPM 中建立了表面張力IPF模型,研究液滴變形和流體在固壁上的潤(rùn)濕行為.Ma 等[75]用“凍結(jié)粒子”簇表達(dá)固體顆粒,應(yīng)用MaPPM 對(duì)包含1024 個(gè)顆粒的氣固懸浮系統(tǒng)進(jìn)行了高分辨率直接數(shù)值模擬,并研究了氣固系統(tǒng)動(dòng)態(tài)多尺度結(jié)構(gòu).Xiong 等[76-77]用MaPPM 模擬研究了有數(shù)十萬(wàn)個(gè)固體顆粒的流體?顆粒流化過(guò)程.

3 SPH 模擬多相流的多種實(shí)現(xiàn)方式

總結(jié)SPH 模擬多相流的各文獻(xiàn)可發(fā)現(xiàn),SPH 模擬多相流的實(shí)現(xiàn)方式大致可分為4 類(lèi):(1)雙流體模型的拉格朗日求解器;(2)多相SPH 方法;(3) SPH 與其他離散方法耦合;(4) SPH 和基于網(wǎng)格的方法耦合.

3.1 TFM 的拉格朗日求解器

類(lèi)似于SPH 在單相流中的應(yīng)用,SPH 方法也可以作為T(mén)FM 的拉格朗日求解器模擬兩相流動(dòng).

TFM[30-32]作為具有不同界面或分散相(氣泡或固體顆粒)兩相流的連續(xù)方法描述.各相均分別由一套質(zhì)量、動(dòng)量和能量守恒方程單獨(dú)描述.由兩相方程中的相間作用項(xiàng)來(lái)耦合相間質(zhì)量、動(dòng)量和能量傳遞.

對(duì)于流固兩相流,各相控制方程表達(dá)為下列方式[30]

流體相

其中θ是各相體積分?jǐn)?shù),β是阻力系數(shù),下標(biāo)f 和s 分別表示液相和固相.兩組方程之間的約束條件是θf(wàn)+θs=1 和0 ≤θf(wàn),θs≤ 1.耦合兩相的相互作用由動(dòng)量方程中的曳力項(xiàng)來(lái)表示.

與SPH 方法求解單相流動(dòng)類(lèi)似,將兩相離散為兩組獨(dú)立SPH 粒子,用于構(gòu)建TFM 的拉格朗日求解器.而兩相SPH 粒子的分布自然滿(mǎn)足約束條件θf(wàn)+θs=1 和0 ≤θf(wàn),θs≤ 1;與其他方法類(lèi)似,解決兩相問(wèn)題的關(guān)鍵是相間相互作用及其離散表示.

Monaghan 和Kocharyan[12]用SPH 方法求解了含塵兩相流動(dòng).動(dòng)量方程中的相間相互作用項(xiàng)被離散化為曳力核函數(shù)計(jì)算

其中η2是防止r2ij=0 時(shí)出現(xiàn)奇點(diǎn)的約為～0.001h2的剪切常數(shù).曳力系數(shù)β與表示空氣和塵埃的SPH粒子屬性相關(guān),并表達(dá)為

其中dp是塵埃直徑,CD是與顆粒雷諾數(shù)Rep=ρgμg|uf?us|/dp有關(guān)的標(biāo)準(zhǔn)阻力系數(shù).兩個(gè)SPH 粒子沿中心線(xiàn)動(dòng)量和角動(dòng)量保持守恒.

Monaghan[78]發(fā)展了一種隱式曳力模型,以提高SPH 計(jì)算的準(zhǔn)確性和穩(wěn)定性.粒子i的速度基于曳力引起的減速度更新

其中上標(biāo)0 表示應(yīng)用曳力項(xiàng)前的初始值,Δt為時(shí)間步長(zhǎng),系數(shù)sij為

Laibe 和Price[79]采用可變光滑長(zhǎng)度的曳力核函數(shù)中來(lái)進(jìn)行氣體?粉塵混合物的模擬.Xiong 等[76]在液固懸浮體系模擬中引入了與固體體積分?jǐn)?shù)相關(guān)的固相黏度系數(shù)計(jì)算式,

3.2 多相SPH 方法

SPH 方法在流動(dòng)模擬中顯示出巨大的潛力,研究人員擴(kuò)展SPH 方法用于多相流模擬,建立多相 SPH方法[15,21-29].在多相SPH 方法中,兩個(gè)不同的相被離散成兩組具有不同屬性的SPH 粒子,并且相間相互作用由不同的SPH 參數(shù)隱式描述.它是 SPH 模擬兩相流的最常用方法.與TFM 的拉格朗日求解器不同,多相SPH 方法采用隱式相互作用來(lái)耦合兩相,并將固體/氣泡視為SPH 粒子的集合.

多相SPH 方法是SPH 方法用于多相流模擬的自然延伸.它包含SPH 方法在單相流模擬中的特點(diǎn)和優(yōu)勢(shì),并引入了各種兩相相互作用(如表面張力)模型來(lái)完成多相流模擬,可精細(xì)準(zhǔn)確地描述界面行為.

3.2.1 液?液體系

液?液體系的模擬是多相SPH 較早期及較常見(jiàn)的應(yīng)用.液?液系統(tǒng)中兩種流體的特性,例如密度和黏度,通常相差不大.

瑞利?泰勒(Rayleigh?Taylor,R-T)不穩(wěn)定性是在具有不同性質(zhì)的兩種液體之間界面上發(fā)展和演變的不穩(wěn)定性.Cummins 和Rudman[56]引入人工表面張力首次用多相SPH 方法進(jìn)行 R-T 不穩(wěn)定性模擬.Tartakovsky 和Meakin[27]發(fā)展了基于數(shù)密度的多相SPH,消除了兩種流體之間的人工表面張力,模擬得到更真實(shí)的界面行為和R-T 不穩(wěn)定性.Hu 和Adams[58]發(fā)展多相投影方法模擬了不可壓不互溶液?液兩相流動(dòng);還發(fā)展基于恒定密度的 SPH 投影方法,并采用CSF 表面張力模型模擬了R-T 不穩(wěn)定性[80].

開(kāi)爾文?亥姆霍茲(Kelvin-Helmholtz,K-H)不穩(wěn)定性是不同水平速度的兩種平行流體間的另一種不穩(wěn)定性現(xiàn)象.Shadloo 和Yildiz[81]應(yīng)用多相SPH方法模擬研究了兩種不可壓不互溶液體的K-H 不穩(wěn)定性.

Shao[82]在界面處應(yīng)用壓力連續(xù)性條件,用多相SPH 方法解耦求解液液兩相.Tong 和Browne[28]采用基于恒定密度的SPH 投影方法研究?jī)煞N流體之間的傳熱和界面處的Marangoni 力效應(yīng).

3.2.2 氣?液體系

相比于液?液體系,氣?液體系中兩相密度比會(huì)高得多,可以達(dá)到數(shù)百或數(shù)千量級(jí).處理高密度比是SPH 方法模擬氣液體系的一個(gè)挑戰(zhàn).常規(guī)SPH 方法還無(wú)法成功處理密度比高于10 的兩相流體[83].SPH 方法進(jìn)行氣液體系模擬的另一個(gè)困難是消除高密度比下表面張力導(dǎo)致的壓力變化.

Colagrossi 和Landrini[83]發(fā)展了多相SPH 方法中兩相界面表征的校正方法,并模擬了帶空氣環(huán)境的潰壩過(guò)程.Das 和Das[84]采用CSF 模型和核函數(shù)近似校正方式發(fā)展多相SPH 方法,成功模擬了水池中氣體從氣孔注入和氣泡演化的過(guò)程.

Szewc 等[26,85]采用密度校正和表面張力CSF模型發(fā)展多相ISPH,模擬了兩相密度比為1000 的液體中氣泡上升變化過(guò)程,并獲得了氣?液體系中氣泡曳力系數(shù)與雷諾數(shù)的關(guān)系.Lind 等[62]采用多相ISPH 方法模擬了密度比為1000 的氣體中液?滴變形過(guò)程.

孫鵬楠等[86]應(yīng)用粒子平移技術(shù)(particle shifting technique,PST)[87-88]改善CSF 模型計(jì)算時(shí)的界面失穩(wěn),成功模擬了氣泡在水中上浮及破碎過(guò)程.

3.2.3 液?固體系

在應(yīng)用多相SPH 方法模擬的液?固體系時(shí),固體常被處理為流動(dòng)的邊界或性質(zhì)截然不同的剛體SPH 粒子.

Akinci 等[89]提出了一種在多相 SPH 方法中通過(guò)流體力雙向耦合流體和剛性物體的方法,實(shí)現(xiàn)了船只通過(guò)橋梁和艦艇在波濤洶涌的大海中航行的模擬.

Tofighi 等[90]將固體表達(dá)為具有極高黏度的SPH 粒子的集合,用于模擬流體中剛體的運(yùn)動(dòng)過(guò)程.Fourtakas 和Rogers[91]提出了一個(gè)液?固體系多相SPH 模型,該模型結(jié)合了液體和固體之間黏附效應(yīng)涉及的屈服變形、剪切和懸浮等作用,用于易蝕壩潰裂過(guò)程的模擬.

3.2.4 氣?固體系

Ma 等[75]采用MaPPM 方法實(shí)現(xiàn)氣固懸浮體系的高分辨率直接數(shù)值模擬.在MaPPM 中,固體顆粒視作一組“凍結(jié)的SPH 粒子”.Xiong 等[76-77]采用MaPPM 對(duì)氣固流態(tài)化進(jìn)行了大規(guī)模直接數(shù)值模擬,以亞網(wǎng)格的精度模擬了多達(dá)30 240 個(gè)固體顆粒和超過(guò)10 億個(gè)流體SPH 粒子,該模擬基于并行GPU 完成.

3.3 耦合SPH 和其他離散方法

在多相流模擬中,SPH 方法通常適合求解連續(xù)相或變形較大的相.SPH 很容易與其他離散方法耦合來(lái)模擬流體?固體系統(tǒng)或流體中的可變形物體.在耦合方法中,SPH 和其他離散方法相結(jié)合,在其特定的適用領(lǐng)域中獲益.這種耦合將SPH 方法的應(yīng)用擴(kuò)展到具有不同特性的多相流;此外,也將不同拉格朗日方法的優(yōu)點(diǎn)結(jié)合到此類(lèi)模擬中.

對(duì)于流體?顆粒體系,應(yīng)用離散元方法(discrete element method,DEM)跟蹤每個(gè)固體顆粒的運(yùn)動(dòng)是很自然的做法.SPH 粒子和DEM 粒子的相互作用(曳力和浮力)表達(dá)為接觸的兩相粒子之間物理量差值的函數(shù).SPH 和DEM 的耦合以亞顆粒的精度實(shí)現(xiàn)了流體?顆粒懸浮體系的直接數(shù)值模擬.

Potapov 等[92]建立SPH-DEM 耦合方法,模擬研究了圓柱繞流和兩平行剪切板之間的流?固體系二維流動(dòng).Huang 和Nydal[93]結(jié)合SPH 和DEM 研究固體顆粒在頂蓋驅(qū)動(dòng)方腔流中的運(yùn)動(dòng)過(guò)程.Robinson等[25]應(yīng)用SPH-DEM 方法模擬了三維液柱中的顆粒沉降,流體對(duì)固體顆粒的作用力表示為

其中Vp是顆粒的體積,fd為流體作用在顆粒上的曳力,其取決于局部的固相體積分?jǐn)?shù)θs和流體?固體顆粒之間的滑移速度uslip=uf?ui.壓力梯度和應(yīng)力張量采用SPH 插值方法計(jì)算

SPH 粒子上的相間力由每個(gè)DEM 粒子貢獻(xiàn)量的加權(quán)平均得到

El Shamy 和Sizkow[94]將SPH 方法和DEM 方法耦合用于評(píng)估動(dòng)態(tài)激勵(lì)下飽和顆粒土壤的液化過(guò)程.Robb 等[95]耦合SPH-DEM 方法模擬了不規(guī)則浮冰接近靜止覆蓋物和上游浮冰在障礙物上積聚等過(guò)程,建立浮冰動(dòng)力學(xué)研究的定量工程工具.Liu 等[96]和Ji 等[97]在基于多面體單元的DEM 模型中引入閔可夫斯基加和操作簡(jiǎn)化不規(guī)則固體顆粒間的碰撞作用計(jì)算,并耦合SPH 方法建立了液固多相流動(dòng)模擬方法,成功模擬了包含石塊運(yùn)動(dòng)的潰壩過(guò)程以及石塊落水過(guò)程.

Polwaththe-Gallage 等[98]使用基于DEM 的彈簧網(wǎng)絡(luò)模型來(lái)表示紅細(xì)胞的黏彈性細(xì)胞膜,并采用SPH 方法對(duì)細(xì)胞內(nèi)外的流體進(jìn)行建模.耦合模型用于模擬毛細(xì)血管中細(xì)胞的變形和運(yùn)動(dòng),并研究了細(xì)胞間相互作用對(duì)細(xì)胞運(yùn)動(dòng)的影響.

最近,陳飛國(guó)和葛蔚[99]將表征氣體運(yùn)動(dòng)的擬顆粒模型(pseudo-particle model,PPM)[100-101]和光滑粒子流體動(dòng)力學(xué)耦合,對(duì)液滴破碎和油膜霧化等氣液流動(dòng)過(guò)程進(jìn)行了高精度模擬.與多相SPH 方法中僅用不同參數(shù)區(qū)分氣液兩相的方式不同,PPM-SPH 耦合方法中由截然不同的模型表達(dá)氣液兩相,更能反映兩相行為差異.

3.4 耦合SPH 和基于網(wǎng)格的方法

SPH 方法在流體流動(dòng)模擬中非常成功,特別是在界面流動(dòng)、大變形流動(dòng)和多相流動(dòng)中.然而,為了獲得高精度模擬結(jié)果,需要采用大量的計(jì)算粒子表示流體,對(duì)計(jì)算要求極高.而對(duì)于單相的流動(dòng)主體,完全的拉格朗日求解并不必需也不高效.基于上述認(rèn)識(shí),將SPH 方法和基于網(wǎng)格的方法耦合可以獲得精度和效率間的平衡.這些耦合方法的關(guān)鍵為重疊區(qū)域中兩類(lèi)方法的銜接處理.

Marrone 等[102]發(fā)展SPH 和有限體積CFD(computational fluid dynamics)的耦合方法,模擬了存在變形和破碎的自由表面流.流動(dòng)主體部分采用結(jié)構(gòu)化網(wǎng)格求解,變形較大的部分采用SPH 方法求解.Napoli 等[103]考慮了CFD 和SPH 重疊區(qū)域的物性匹配,提高了SPH-CFD 耦合方法模擬的準(zhǔn)確性.Deng 等[104]將SPH 和有限體積CFD 耦合用于求解氣固流動(dòng)雙流體模型.分散的固體顆粒進(jìn)行類(lèi)似流體的連續(xù)處理,并由基于GPU 的SPH 計(jì)算求解;而流體運(yùn)動(dòng)模擬則采用有限體積CFD 在Ansys FLUENT軟件中實(shí)施.

Feng 等[105]耦合SPH 和有限元方法(finite element method,FEM)模擬了磨料水射流切割的單次加速過(guò)程.Zhang 等[106]在 SPH 中加入了大渦模擬(large eddy simulation,LES)表征湍流效應(yīng),并使用有限元虛擬邊界法描述固體顆粒,模擬了流體中顆粒沉降過(guò)程,驗(yàn)證了SPH 改進(jìn)方法用于多相流模擬研究的能力.

SPH-FEM 耦合模型常用于研究流?固耦合問(wèn)題,SPH 的引入降低了流體?固壁間耦合處理復(fù)雜性[107],而FEM 可以提高復(fù)雜固壁的處理能力以及計(jì)算準(zhǔn)確性和穩(wěn)定性[108-109].Long 等[110]在SPHFEM 耦合模型中加入共軛傳熱計(jì)算,用于解決熱流構(gòu)耦合問(wèn)題.

4 多相流 SPH 方法的其他討論

4.1 近似修正方法

4.1.1 核近似修正

Chen 等[111]提出了一種修正光滑粒子方法(corrective smoothed-particle method,CSPM),通過(guò)提高核近似的精度來(lái)解決張力不穩(wěn)定性問(wèn)題

CSPM 已成功用于固壁邊界[111]和兩相界面處[28,84]的流體改進(jìn)求解.

4.1.2 粒子近似修正

在多相體系中直接應(yīng)用質(zhì)量密度表達(dá)式可能會(huì)引起人工表面張力并扭曲界面處的密度計(jì)算[83].Tartakovsky 和Meakin[112]提出使用粒子數(shù)密度來(lái)消除不同相屬性對(duì)密度計(jì)算的影響,從而避免人工表面張力.粒子數(shù)密度計(jì)算公式為

此外,這種粒子數(shù)密度近似方法適用于高密度比的多相流動(dòng)模擬.

4.1.3 一致性修正

SPH 近似表達(dá)中的一致性(Cn一致性)表示為可由SPH 方法精確再現(xiàn)多項(xiàng)式的階數(shù)[113].它受光滑核函數(shù)、光滑長(zhǎng)度、邊界條件和粒子分布等因素影響.通常,傳統(tǒng)SPH 方法可認(rèn)為恢復(fù)到一階一致性(C1一致性).但是對(duì)于具有不均勻粒子分布、光滑長(zhǎng)度與顆粒間距不同或邊界(固壁邊界或相界面)處,傳統(tǒng)SPH 方法甚至不能完全滿(mǎn)足零階(C0)一致性[114].

Liu 和Liu[36,114] 提出了有限粒子法(finite particle method,FPM)來(lái)提高SPH 方法的一致性.在FPM 中,場(chǎng)函數(shù)f(x)及其導(dǎo)數(shù)的近似值修正為[113]

對(duì)于遇到邊界缺陷和粒子分布敏感性問(wèn)題的情況,FPM 大大提高了SPH 方法的一致性.但對(duì)于粒子分布高度無(wú)序的情況,由于病態(tài)校正矩陣導(dǎo)致數(shù)值求解不穩(wěn)定,FPM 模擬可能也會(huì)崩潰[115].合理假設(shè)粒子自我導(dǎo)向的貢獻(xiàn)在光滑核函數(shù)的導(dǎo)數(shù)中占主導(dǎo)地位,Zhang 等[106]將修正矩陣簡(jiǎn)化為如下對(duì)角矩陣建立解耦有限粒子法(decoupled finite particle method,DFPM)

因此,f(x)及其導(dǎo)數(shù)的粒子近似為

4.2 張力不穩(wěn)定性

張力不穩(wěn)定性是應(yīng)用SPH 方法時(shí)一個(gè)眾所周知的問(wèn)題,它會(huì)導(dǎo)致粒子過(guò)密聚集或粒子被炸開(kāi)[33].張力不穩(wěn)定性在SPH 模擬中時(shí)常被觀(guān)察到,一直是早期SPH 研究[116]的關(guān)注焦點(diǎn).即使壓力為正值,張力不穩(wěn)定性仍會(huì)發(fā)生.Swegle 等[117]和Monaghan[118]給出了應(yīng)力狀態(tài)穩(wěn)定性或不穩(wěn)定性的判據(jù);應(yīng)力不穩(wěn)定增長(zhǎng)的條件為W′′

這里 是光滑核函數(shù)的二階導(dǎo)數(shù),S是應(yīng)力狀態(tài).該條件表明,核函數(shù)對(duì)SPH 方法的穩(wěn)定模擬至關(guān)重要.核函數(shù)的二階導(dǎo)數(shù)在整個(gè)支持域中應(yīng)該滿(mǎn)足非負(fù)性.同時(shí),建議采用雙曲型核函數(shù)來(lái)消除張力不穩(wěn)定性[39].

4.3 高雷諾數(shù)流動(dòng)

眾所周知,與其他粒子方法類(lèi)似,傳統(tǒng)SPH 方法不適合模擬高雷諾數(shù)下的流動(dòng),這是由于高渦度條件下會(huì)發(fā)生張力不穩(wěn)定性以及在負(fù)壓條件下會(huì)產(chǎn)生非真實(shí)的粒子空缺[119].

Sun 等[119-120]在δ-SPH 方法[121]中實(shí)施張力不穩(wěn)定性控制(tensile instability control,TIC)和粒子平移技術(shù)[87-88]建立δ+-SPH 方法.在δ+-SPH 方法中,連續(xù)性方程引入了一個(gè)額外的數(shù)值擴(kuò)散項(xiàng),并將粒子平移重構(gòu)粒子分布.δ+-SPH 方法克服了SPH 方法難以模擬高雷諾數(shù)流動(dòng)的局限性.

4.4 表面張力系數(shù)

如前面小節(jié)里所述,SPH 方法中兩相界面處的表面張力主要由CSF,IPF 或CPM 等模型表達(dá).CSF 模型直接用實(shí)際表面張力系數(shù)來(lái)表達(dá)宏觀(guān)表面張力.然而,對(duì)于IPF 和CPM 模型,表面張力系數(shù)隱含在模型參數(shù)中,其理論關(guān)系式尚未獲得[27].

Zhou 等[29]采用如下IPF 模型表征兩不互溶液體間的表面張力作用

研究了不同參數(shù)s下固體基質(zhì)上液滴在液體中的變形情況.兩相表面張力系數(shù)γ則可由躺滴法根據(jù)下式計(jì)算得到

其中Δρ是兩種液體之間的密度差,H是液滴的半高(見(jiàn)圖7).模型系數(shù)s和物理表面張力系數(shù)γ之間的關(guān)系表達(dá)如下[29]

圖7 躺滴示意Fig.7 Profile of a sessile drop

Xu 等[122]和Yang 等[33]應(yīng)用CPM 模型模擬了不同半徑的黏性液滴的粒子分布情況,而表面張力系數(shù)γ可基于拉普拉斯方程計(jì)算得到

其中pl是液滴中心壓力,pg是遠(yuǎn)離液滴的蒸汽壓,R是液滴的半徑.液滴中心壓力則根據(jù)液滴中心粒子密度由范德華狀態(tài)方程計(jì)算得到,不同半徑的液滴中心壓力與1/R線(xiàn)性擬合得到的斜率即為表面張力系數(shù),如圖8 所示.

圖8 不同液滴半徑下的液滴中心壓力Fig.8 Pressure at the center of drops for various drop radii

5 結(jié)論

作為無(wú)網(wǎng)格方法的代表,SPH 在過(guò)去的幾十年中發(fā)展迅速,并在許多研究領(lǐng)域得到了廣泛的應(yīng)用.由于其天然的拉格朗日特征,SPH 方法具有處理界面、大變形、不連續(xù)性和多物理場(chǎng)的能力,在多相流的模擬研究中顯示出巨大的潛力.許多研究者在SPH 方法模擬多相流方面已獲得眾多成果.本工作試圖從實(shí)現(xiàn)方式的異同歸類(lèi)這些多相流SPH 模擬研究.

(1)鑒于SPH 方法成功應(yīng)用于單相流動(dòng)模擬,作為T(mén)FM 的完全拉格朗日求解器則是SPH 方法模擬多相流的直接方法.該方式需要多相系統(tǒng)各相耦合項(xiàng)以及建模參數(shù)(例如壓力和黏度)的顯式表達(dá)式.

(2)多相SPH 方法則對(duì)每個(gè)相分別SPH 建模,并且自然地將具有不同參數(shù)的兩個(gè)相耦合.多相SPH 方法是SPH 方法中應(yīng)用最多的多相流模擬方法,已成功應(yīng)用于各種不同體系的研究.

(3) SPH 方法還能與其他離散方法(如分子動(dòng)力學(xué)、耗散粒子動(dòng)力學(xué)和離散單元法)相耦合模擬各類(lèi)多相流.耦合方法保持了每相對(duì)應(yīng)的不同離散方法的特征和優(yōu)勢(shì).

(4)此外,為了降低使用大量SPH 粒子進(jìn)行多相流模擬的計(jì)算成本,還可以使用基于網(wǎng)格的方法求解較為簡(jiǎn)單的主體流動(dòng),而僅在界面處應(yīng)用SPH 方法.歐拉?拉格朗日耦合方法為多相流模擬提供了準(zhǔn)確性和效率之間的折衷方案,并顯示出模擬實(shí)際多相流體系的巨大潛力.

然而,為了提高多相流 SPH 模擬的準(zhǔn)確性和效率,仍有一些問(wèn)題需要解決.例如,(1) SPH 參數(shù)與流體物性之間的關(guān)系.不同的加權(quán)核函數(shù)會(huì)導(dǎo)致不同的有效物性[123-124].為了穩(wěn)定模擬,還引入了如人工黏度[125]等人工參數(shù).因此,模型參數(shù)和實(shí)際物性的校準(zhǔn)會(huì)隨著不同核函數(shù)、光滑長(zhǎng)度和近似方式而有所差異.(2)計(jì)算成本高的解決方案.SPH 方法雖然適用于模擬具有大變形和分明相界面的流動(dòng),但其計(jì)算成本非常高,尤其是模擬多相流.一種解決方案是實(shí)現(xiàn)大規(guī)模并行計(jì)算[22,77,91,126-127] 以及基于GPU 等加速計(jì)算技術(shù)[128-134];另一種方案是減少不必要的計(jì)算,如耦合其他離散方法表達(dá)分散體,耦合基于網(wǎng)格方法模擬連續(xù)相主體平緩流動(dòng)等.(3)實(shí)驗(yàn)驗(yàn)證和增強(qiáng)模擬.有許多實(shí)驗(yàn)驗(yàn)證了SPH 模擬,例如潰壩[64,135]和氣?液霧化[22].此外,SPH 方法還廣泛應(yīng)用于計(jì)算機(jī)視覺(jué)用于顯示真實(shí)的流體行為[60,136].實(shí)驗(yàn)、模擬和視覺(jué)技術(shù)的結(jié)合可能會(huì)激勵(lì)在計(jì)算機(jī)上進(jìn)行現(xiàn)實(shí)中困難或無(wú)法實(shí)現(xiàn)的增強(qiáng)虛擬實(shí)驗(yàn).