陳飛國 * 葛 蔚

* (中國科學院過程工程研究所多相復(fù)雜系統(tǒng)國家重點實驗室,北京 100190)

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

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

引言

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

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

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

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

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

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

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

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

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

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

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

1 SPH 方法

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

1.1 SPH 方法的基本原理

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

1.1.1 核近似

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

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

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

函數(shù)導數(shù)由核函數(shù)近似表達為

1.1.2 粒子近似

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

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

其在粒子i處的值為

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

1.2 控制方程的離散化

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

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

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

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

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

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

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

1.3 狀態(tài)方程

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

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

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

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

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

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

1.4 核函數(shù)

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

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

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

(2)緊致性

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

(3)非負性

(4)當光滑長度趨于零時,光滑核函數(shù)滿足狄拉克δ函數(shù)(δ函數(shù)性)

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

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

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

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

1.5 表面張力

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

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

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

1.5.1 CSF 模型

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

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

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

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

1.5.2 IPF 模型

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

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

文獻中常用的IPF 模型包括Lennard?Jones 力模型[46]

更多的IPF 模型見于Wang 等[15]的綜述.

1.5.3 CPM

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

1.6 固體邊壁

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

1.6.1 粒子表示法

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

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

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

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

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

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

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

1.6.2 幾何表示法

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

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

2 多相流模擬其他SPH 類方法

2.1 不可壓縮SPH

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

中間速度場u*是在沒有壓力梯度項的情況下在時間上向前積分動量方程計算得到.壓力泊松方程基于中間速度場獲得并求解獲得達到不可壓縮條件的壓力場

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

2.2 移動粒子半隱式法

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

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

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

2.3 宏觀擬顆粒方法

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

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

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

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

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

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

3.1 TFM 的拉格朗日求解器

類似于SPH 在單相流中的應(yīng)用,SPH 方法也可以作為TFM 的拉格朗日求解器模擬兩相流動.

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

對于流固兩相流,各相控制方程表達為下列方式[30]

流體相

其中θ是各相體積分數(shù),β是阻力系數(shù),下標f 和s 分別表示液相和固相.兩組方程之間的約束條件是θf+θs=1 和0 ≤θf,θs≤ 1.耦合兩相的相互作用由動量方程中的曳力項來表示.

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

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

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

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

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

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

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

3.2 多相SPH 方法

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

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

3.2.1 液?液體系

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

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

開爾文?亥姆霍茲(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 投影方法研究兩種流體之間的傳熱和界面處的Marangoni 力效應(yīng).

3.2.2 氣?液體系

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

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

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

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

3.2.3 液?固體系

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

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

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

3.2.4 氣?固體系

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

3.3 耦合SPH 和其他離散方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1 近似修正方法

4.1.1 核近似修正

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

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

4.1.2 粒子近似修正

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

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

4.1.3 一致性修正

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

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

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

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

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

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

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

4.3 高雷諾數(shù)流動

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

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

4.4 表面張力系數(shù)

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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