多相流動的光滑粒子流體動力學方法研究綜述1)

陳飛國 * 葛 蔚

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

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

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

引言

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

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

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

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

2)自發捕獲和跟蹤間斷面,而這是使用基于網格的方法難以實現的;

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

4) SPH 粒子大小可靈活調節,對應不同的數值精度要求;

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

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

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

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

1 SPH 方法

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

1.1 SPH 方法的基本原理

簡而言之,SPH 的基本思想[8]可概括如下:用一組帶有物理性質的粒子來離散化問題域,并將流體的控制方程轉變為具有各種邊界條件的常微分方程;然后,采用核函數插值近似的權重積分獲得數值解.SPH 方法的關鍵是插值近似,包括核近似和粒子近似.

1.1.1 核近似

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

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

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

函數導數由核函數近似表達為

1.1.2 粒子近似

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

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

其在粒子i處的值為

這樣,函數的取值及其導數值可以通過離散粒子的加權平均得到.

1.2 控制方程的離散化

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

其中u為速度,ρ為密度,S為應力張量,g為單位質量體積力,e為比內能,q為熱通量.應力張量S采用顯式形式

這里I表示單位張量,p為內壓,剪切應力σ表達為

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

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

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

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

1.3 狀態方程

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

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

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

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

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

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

1.4 核函數

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

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

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

(2)緊致性

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

(3)非負性

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

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

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

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

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

1.5 表面張力

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

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

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

1.5.1 CSF 模型

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

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

通常由下式定義的色函數區分不同的相[37,41]

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

1.5.2 IPF 模型

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

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

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

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

1.5.3 CPM

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

1.6 固體邊壁

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

1.6.1 粒子表示法

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

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

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

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

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

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

圖5 動態粒子模型Fig.5 Schematic illustration of dynamics particles

1.6.2 幾何表示法

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

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

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

2.1 不可壓縮SPH

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

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

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

2.2 移動粒子半隱式法

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

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

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

2.3 宏觀擬顆粒方法

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

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

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

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

3 SPH 模擬多相流的多種實現方式

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

3.1 TFM 的拉格朗日求解器

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

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

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

流體相

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

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

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

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

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

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

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

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

3.2 多相SPH 方法

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

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

3.2.1 液?液體系

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

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

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

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

3.2.2 氣?液體系

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

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

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

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

3.2.3 液?固體系

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

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

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

3.2.4 氣?固體系

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

3.3 耦合SPH 和其他離散方法

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

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

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

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

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

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

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

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

3.4 耦合SPH 和基于網格的方法

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

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

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

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

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

4.1 近似修正方法

4.1.1 核近似修正

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

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

4.1.2 粒子近似修正

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

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

4.1.3 一致性修正

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

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

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

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

4.2 張力不穩定性

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

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

4.3 高雷諾數流動

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

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

4.4 表面張力系數

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

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

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

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

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

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

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

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

5 結論

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

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

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

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

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

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