修正壓力與表面張力計(jì)算的兩相流自由界面運(yùn)動(dòng)模擬

朱曉臨，周韻若，何紅虹

修正壓力與表面張力計(jì)算的兩相流自由界面運(yùn)動(dòng)模擬

朱曉臨，周韻若，何紅虹

(合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，安徽 合肥 230009)

在利用光滑粒子流體幼力學(xué)(SPH)方法進(jìn)行多相流的模擬中，多相流在交界面處由于密度不連續(xù)、粒子界面的壓力計(jì)算出現(xiàn)誤差，出現(xiàn)界面處壓力振蕩，界面破碎的問(wèn)題。針對(duì)上述問(wèn)題，提出了新的壓力梯度近似公式和改進(jìn)的界面處人工斥力公式，使交界面更加清晰、光滑、無(wú)穿透、模擬效果更好。此外，通過(guò)給出基于密度權(quán)重的色函數(shù)計(jì)算公式，對(duì)大密度比多相流界面處的表面張力計(jì)算公式進(jìn)行了改進(jìn)，使得多相流交界面密度過(guò)渡更加光滑，模擬效果更好。最后，通過(guò)潰壩、Rayleigh-Taylor界面不穩(wěn)定性，非Boussinesq鎖定交換問(wèn)題3個(gè)仿真實(shí)驗(yàn)，得到了不同時(shí)刻界面粒子分布圖、界面鋒面距離等。結(jié)果驗(yàn)證了新的壓力梯度近似公式以及界面處人工斥力公式的合理性；通過(guò)空氣中液滴形成仿真實(shí)驗(yàn)，得到了圓形液滴形成的粒子變化圖，結(jié)果驗(yàn)證了該改進(jìn)的表面張力計(jì)算方法的有效性。

SPH方法；多相流；壓力梯度近似公式；人工斥力；色函數(shù)；表面張力

2005年，柳有權(quán)等[1]提出流動(dòng)系統(tǒng)中存在多種不互溶的流體，稱(chēng)為混合流；而若相態(tài)多于一個(gè)，且各相態(tài)之間不互相滲透或發(fā)生化學(xué)反應(yīng)，則稱(chēng)為多相流。近年來(lái)，多相流動(dòng)的數(shù)值模擬技術(shù)是計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics，CFD)領(lǐng)域研究的熱點(diǎn)和前沿課題之一。液液兩相流動(dòng)與氣液兩相流動(dòng)在日常生活中十分普遍，如通過(guò)預(yù)先模擬實(shí)現(xiàn)藥物在人體血管中的精確定位打靶技術(shù)以及氣泡上升至自由表面最終破裂的現(xiàn)象。

而在兩相流動(dòng)數(shù)值模擬技術(shù)中，如何處理兩相交界面處的界面穩(wěn)定性是最近幾年研究中的熱點(diǎn)問(wèn)題。

常見(jiàn)的界面不穩(wěn)定性有2類(lèi)：① 2種不同密度的流體在重力或慣性力的作用下的界面不穩(wěn)定性問(wèn)題，通常稱(chēng)之為Rayleigh-Talor界面不穩(wěn)定性問(wèn)題；② 2種流體之間有切向速度差，因此造成界面上的不穩(wěn)定性，通常稱(chēng)之為Kelvin- Helomholtz不穩(wěn)定性。

2003年，COLAGROSSI和LANDRINI[2]將1995年MONAGHAN和KOCHARYAN[3]提出的模擬小密度比兩相流的方法應(yīng)用于氣液兩相流模擬，發(fā)現(xiàn)交界面處出現(xiàn)劇烈的數(shù)值振蕩。該方法使用弱可壓縮模型建模，根據(jù)粒子的體積而不是密度進(jìn)行建模，進(jìn)行壓力梯度計(jì)算，避免了氣-水界面密度不連續(xù)的問(wèn)題。此方法需要極小的時(shí)間步長(zhǎng)作為支撐，整個(gè)模擬需要的時(shí)間較長(zhǎng)。

2006年，HU和ADAMS[4]提出了可計(jì)算宏觀(guān)和微觀(guān)流場(chǎng)的光滑流體粒子動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)方法，采用不同的粒子近似函數(shù)，即粒子近似函數(shù)是針對(duì)鄰近點(diǎn)體積而不是密度，但其并未計(jì)算大密度比復(fù)雜界面流場(chǎng)。

2007年，HU和ADAMS[5]提出了一種不可壓縮的多相SPH方法，主要通過(guò)一種新的多相投影公式來(lái)離散梯度與散度算子，但是該方法需要求解壓力泊松方程，整個(gè)計(jì)算過(guò)程十分耗時(shí)。

2008年，SOLENTHALER和PAJAROLA[6]提出一種使用標(biāo)準(zhǔn)SPH模型去模擬大密度比兩相流。其主要通過(guò)計(jì)算粒子數(shù)密度而不是計(jì)算粒子密度去保證兩相界面密度的連續(xù)性，從而推出壓力項(xiàng)與黏性項(xiàng)的計(jì)算公式。該方法簡(jiǎn)單，但仍存在局限性。其直接基于密度進(jìn)行修正的方法有效解決了兩相流密度上的不連續(xù)性，但是在計(jì)算壓力項(xiàng)，尤其針對(duì)體積較小的流體時(shí)，兩相流體間還是會(huì)存在較大的間隙，無(wú)法保證交界面的連續(xù)性。

2013年，MONAGHAN和RAFIEE[7]通過(guò)在交界面上引入斥力項(xiàng)維持界面的光滑與完整。該方法簡(jiǎn)單，但是對(duì)于聲速的使用依舊是非物理的，同時(shí)將界面處粒子視為剛性邊界粒子，對(duì)粒子的運(yùn)動(dòng)范圍設(shè)置限制，其好處是穩(wěn)定界面，但是卻是不符合物理規(guī)律，理論上界面上的粒子可以向各個(gè)方向進(jìn)行運(yùn)動(dòng)，模擬缺乏處理湍流的能力。

2015年，CHEN等[8]針對(duì)密度差較大的多相流，提出了自己的SPH模型，其舍棄了文獻(xiàn)[4]提出的顆粒間壓力的想法，通過(guò)將界面粒子的鄰近粒子中非同相粒子看成同相粒子，其速度、體積、位置等信息與初始異相粒子相同。質(zhì)量與密度使用當(dāng)前相流體粒子的質(zhì)量與密度；同時(shí)提出了一種新的針對(duì)大密度比的多相流密度初始化的方法，以及為避免初始負(fù)壓強(qiáng)，在狀態(tài)方程中加入初始?jí)簭?qiáng)項(xiàng)；但是該方法并未在界面處加入人工斥力項(xiàng)，界面較為粗糙。

2015年，LIND等[9]針對(duì)氣-水界面的大密度比問(wèn)題提出一種不可壓縮-可壓縮建模方法。將水相作為不可壓縮相，空氣相作為弱可壓縮相進(jìn)行建模，使得模擬中保證了界面處材料的不連續(xù)性以及壓力和運(yùn)動(dòng)速度的連續(xù)性。但是該方法在求解水相的壓力時(shí)，同樣需要解費(fèi)時(shí)的壓力泊松方程，計(jì)算量較大。

2018年，KRIMI等[10]提出加入連續(xù)壓力表面方程的SPH模型模擬界面多相流，主要加入一個(gè)壓力表面系數(shù)，提高了界面處的穩(wěn)定性與光滑性；但是該方法在模擬過(guò)程中，對(duì)于界面系數(shù)的選擇較為困難，缺乏量化的標(biāo)準(zhǔn)去選擇表面系數(shù)。

2019年，楊秋足等[11]提出了一種利用弱可壓縮SPH方法來(lái)對(duì)大密度比多相流模擬建模的方法；主要在求解界面壓力與密度處引入了黎曼解的相關(guān)思想，但是此方法計(jì)算復(fù)雜，耗時(shí)較長(zhǎng)。

盡管經(jīng)典的SPH方法成功地模擬了單相流體的流動(dòng)，但由于多相流體的流動(dòng)中多相物理?xiàng)l件的不同導(dǎo)致單相流體的SPH方法難以適用于多相流。主要的問(wèn)題是兩相流密度在界面處是跳躍的，但是SPH方法是一個(gè)平滑過(guò)渡的過(guò)程，每個(gè)粒子都有一個(gè)緊致的支持域，通過(guò)緊致的支持域內(nèi)均勻分布的粒子計(jì)算出當(dāng)前粒子的各項(xiàng)物理屬性。當(dāng)支持域與界面相交時(shí)，界面另一側(cè)的密度會(huì)錯(cuò)誤地影響局部密度和壓力場(chǎng)的計(jì)算，從而影響有關(guān)粒子的加速度與速度的計(jì)算，導(dǎo)致粒子下一個(gè)時(shí)間步長(zhǎng)的位置錯(cuò)誤。因此，本文針對(duì)多相流界面處的密度、壓力、表面張力等物理量進(jìn)行單獨(dú)處理，以保證密度梯度的連續(xù)和表面張力的連續(xù)，避免界面處的粒子在運(yùn)動(dòng)過(guò)程中出現(xiàn)振蕩現(xiàn)象。

1 控制方程與SPH方法

1.1 控制方程

1977年，LUCY[12]首次提出SPH方法，之后被GINGOLD和MONAGHAN[13]用于模擬天體物理問(wèn)題。1991年，LIBERSKY和PETSCHEK[14]將其擴(kuò)展到流體模擬領(lǐng)域，促進(jìn)了SPH對(duì)于水動(dòng)力學(xué)的應(yīng)用研究。

流體運(yùn)動(dòng)控制方程為

利用弱可壓縮模型描述不可壓縮流體，采用如下?tīng)顟B(tài)方程[8]

1.2 SPH方法

其中，為平滑核函數(shù)(smoothing kernel function)，為核函數(shù)的支持域長(zhǎng)度。

其中，為粒子的鄰近粒子；m為粒子的質(zhì)量；為粒子的密度。

分別對(duì)()進(jìn)行梯度和拉普拉斯運(yùn)算，即分別對(duì)其光滑核函數(shù)進(jìn)行梯度和拉普拉斯運(yùn)算

1.3 黏性項(xiàng)

考慮到流體具有黏性，動(dòng)量方程中黏性項(xiàng)的形式為

本文采取文獻(xiàn)[16]提出的黏性項(xiàng)表達(dá)式，即

其中，=-；=-；=；?W=?(-,)|=r。

當(dāng)使用SPH方法進(jìn)行流體運(yùn)動(dòng)的計(jì)算過(guò)程中，需要采用人工黏性項(xiàng)以避免激波附近出現(xiàn)非物理震蕩。本文采用MONAGHAN和GINGOLD[17]提出的人工黏性項(xiàng)，即

2 本文工作

2.1 新的粒子壓力梯度近似公式

針對(duì)上述多相流界面處存在的問(wèn)題，下面對(duì)多相流界面處的壓力梯度近似公式進(jìn)行改進(jìn)，解決多相流數(shù)值模擬中界面不光滑的問(wèn)題，并與文獻(xiàn)[8]的結(jié)果進(jìn)行對(duì)比。文獻(xiàn)[4]提出的界面處的壓力項(xiàng)為

其中，粒子間壓力項(xiàng)為

式(18)對(duì)于小密度比兩相流成立，但是對(duì)于大密度比兩相流，在界面處，p?p，<<，此時(shí)，式(18)變?yōu)?/p>

可以發(fā)現(xiàn)，加入的粒子間壓力項(xiàng)主要是由密度較低的流體的壓力引導(dǎo)。此時(shí)可以等價(jià)為只有小密度流體對(duì)大密度流體有作用，而大密度流體對(duì)小密度流體無(wú)作用，顯然有問(wèn)題。同時(shí)，在界面處，進(jìn)行鄰近粒子搜索時(shí)，混合著兩相粒子，會(huì)對(duì)顆粒間壓力項(xiàng)計(jì)算產(chǎn)生誤差。

文獻(xiàn)[8]對(duì)式(17)進(jìn)行了改進(jìn)，對(duì)界面處的壓力梯度近似公式進(jìn)行了修改，具體形式為

由式(20)可以發(fā)現(xiàn)，在計(jì)算界面處粒子間壓力項(xiàng)時(shí)，使用的均為粒子b的相關(guān)信息。文獻(xiàn)[8]所做的針對(duì)性處理是將粒子a的鄰近粒子中的非同相粒子均人為地看作同相粒子，使用原有的粒子的信息(除了質(zhì)量、密度)，如圖1所示。

上述做法需要在計(jì)算中人為地更改粒子的物理屬性，增加了額外的計(jì)算量。

本文基于牛頓第三定律，作用力等于反作用力，在界面處大密度粒子與小密度粒子對(duì)彼此的作用力應(yīng)是相等的，這樣才可以保證界面的光滑與完整；否則，粒子會(huì)因?yàn)閴毫Σ坏仍诮缑嫣幊霈F(xiàn)振蕩，將式(20)改成對(duì)稱(chēng)形式，保證界面處壓力的連續(xù)，得到新的粒子壓力梯度近似式為

使用當(dāng)前粒子與鄰近粒子的物理信息，避免了人為修改粒子的物理信息，減小了計(jì)算量。對(duì)比式(17)，式(21)解決了大密度粒子對(duì)小密度粒子無(wú)影響的問(wèn)題，同時(shí)保證了界面處壓力的連續(xù)性，其作為粒子壓力梯度近似公式是合理的。

在壓力場(chǎng)中

2.2 改進(jìn)的人工斥力公式

在實(shí)驗(yàn)中，發(fā)現(xiàn)仍存在少量高密度粒子滯留在低密度粒子上方。本文提出了一種新的壓力梯度近似公式，保證低密度粒子和高密度粒子受到的力大小相等，方向相反。但低密度粒子與高密度粒子依舊存在著質(zhì)量上的不連續(xù)，即相互作用力產(chǎn)生的加速度偏差較大，導(dǎo)致出現(xiàn)上述的實(shí)驗(yàn)現(xiàn)象。為了與物理上相斥流體界面現(xiàn)象相對(duì)應(yīng)(即相斥流體界面間存在明顯的分層現(xiàn)象)，故在界面處加入人工斥力，即

式(25)有效解決了不同相粒子在界面處出現(xiàn)混合的情況，保持了界面的光滑性。這與現(xiàn)實(shí)生活中2種密度不同，互不相溶的液體相互接觸的情形相同。

2.3 改進(jìn)的表面張力計(jì)算公式

2010年，ADAMI等[18]提出了多相流中的表面張力計(jì)算公式。

定義色函數(shù)如下

而在實(shí)驗(yàn)中發(fā)現(xiàn)，針對(duì)大密度比多相流，上述的色函數(shù)計(jì)算方式往往會(huì)在界面處出現(xiàn)顏色梯度的跳躍，導(dǎo)致界面出現(xiàn)毛糙等現(xiàn)象。本文給出基于密度權(quán)重的色函數(shù)計(jì)算公式

通過(guò)外法向可以計(jì)算兩相流界面曲率，即

其中，為計(jì)算域的維度。將曲率和外法向代入表面張力公式，得到基于密度權(quán)重的色函數(shù)的改進(jìn)的表面張力計(jì)算式為

根據(jù)上述討論，本文動(dòng)量方程的離散形式為

3 實(shí)驗(yàn)結(jié)果與分析

本節(jié)的前3個(gè)仿真實(shí)驗(yàn)，包括潰壩、R-T不穩(wěn)定性、非Boussinesq交換，檢驗(yàn)本文新的壓力梯度近似公式和改進(jìn)的界面處人工斥力公式的合理性；其中，算例1的密度比大于800，算例2和3的密度比小于2。第4個(gè)仿真實(shí)驗(yàn)——圓形液滴形成，驗(yàn)證本文改進(jìn)的表面張力計(jì)算公式的有效性。

本文的實(shí)驗(yàn)在Windows 10系統(tǒng)上，開(kāi)發(fā)環(huán)境為Visual Studio 2013與openGL，實(shí)驗(yàn)配置為Intel(R) Core(TM) i7-8570H，2.20 GHz CPU，8 G內(nèi)存，NVIDIA Geforce GT 1060顯卡。

3.1 潰壩(算例1)

圖2 二維潰壩模擬粒子效果圖

從圖2可以看出，受重力的作用，初始靜止的水沿著左邊的固壁邊界向右加速運(yùn)動(dòng)(圖2(a1))，水面的高度自左至右不斷降低，水的速度不斷增大，直至到達(dá)右邊的固壁邊界處，與右邊的固壁邊界發(fā)生撞擊，速度減小但壓力增大，水自由面沿墻豎直向上爬升(圖2(a2))。與當(dāng)水到達(dá)最高點(diǎn)且速度降為零時(shí)，重力又使水沿右邊的固壁邊界反向加速運(yùn)動(dòng)，上行與下行的水流在右邊界1/4處匯集，形成了“鼓包”區(qū)域(圖2(a3))，與圖2(b4)相比，本文方法在鼓包形成區(qū)域形態(tài)更為完整，且液體沒(méi)有明顯的壓力振蕩。由于上部的水繼續(xù)下落，受擠壓的鼓包形成了“水射流”(圖2(a4))；“水射流”包裹部分空氣，重力導(dǎo)致水射流向左下方下落，形成類(lèi)似“水射流”沖擊自由表面，撞擊后再次形成“水射流”(圖2(a5))。與圖2(b4)、圖2(b5)相比，本文的方法在界面處粒子無(wú)明顯的“黏滯”狀態(tài)且界面處無(wú)明顯的壓力振蕩。

圖3 三維潰壩模擬粒子效果圖

圖3為三維潰壩實(shí)驗(yàn)結(jié)果。實(shí)驗(yàn)結(jié)果表明，本文方法與文獻(xiàn)[8]的方法結(jié)果進(jìn)行對(duì)比，本文方法在潰壩計(jì)算中的氣-水界面處更加光滑。在流體運(yùn)動(dòng)的過(guò)程中，雖然流體的壓力連續(xù)變化，但是計(jì)算過(guò)程中液面并未出現(xiàn)明顯的壓力振蕩。

圖4給出了液面最高點(diǎn)隨時(shí)間變化的趨勢(shì)。從中可以看出，本文實(shí)驗(yàn)結(jié)果與文獻(xiàn)[18]基本一致，說(shuō)明本文方法在計(jì)算大密度比界面問(wèn)題中的合理性。

3.2 Rayleigh-Taylor界面不穩(wěn)定性(算例2)

采用多相流SPH方法，文獻(xiàn)[5]、文獻(xiàn)[7]和文獻(xiàn)[8]均計(jì)算了R-T界面，其中文獻(xiàn)[8]計(jì)算R-T界面的實(shí)驗(yàn)效果最好。

初始狀態(tài)如圖5(a)所示。當(dāng)=0時(shí)，密度為1 750 kg/m3的流體(紅色)置于密度1 000 kg/m3的流體(藍(lán)色)上方，界面形狀為=1.0-0.15(2p)。取0.05，D=0.003 s。

圖4 水前端與時(shí)間的關(guān)系

圖5 R-T實(shí)驗(yàn)本文方法的不同時(shí)刻粒子分布圖

從本文的實(shí)驗(yàn)結(jié)果(圖5)可以看出，本文算法清晰地捕捉到不同時(shí)刻2種流體界面，界面并未發(fā)現(xiàn)明顯的壓力振蕩，高密度與低密度流體分別呈斜射入狀進(jìn)入對(duì)方區(qū)域。通過(guò)對(duì)最終效果的比較(圖5(d)與圖6)，本文方法在界面處的光滑性更好，計(jì)算中未出現(xiàn)明顯的粒子“集聚”現(xiàn)象，避免了形成虛假的壓力場(chǎng)(即人工表面張力)導(dǎo)致界面曲率變大，界面之間出現(xiàn)狹小的空隙。

圖6 文獻(xiàn)[8]方法

定量驗(yàn)證計(jì)算結(jié)果，圖7給出了低密度流體界面最高點(diǎn)位移隨時(shí)間變化的曲線(xiàn)，從圖6中可以看出，本文實(shí)驗(yàn)結(jié)果基本與LAWRIE 和DALZIEL[19]得到的Layzer理論值重合，同時(shí)與文獻(xiàn)[8]的結(jié)果相比，在范圍內(nèi)，本文得到的-曲線(xiàn)與Layzer理論值更為符合。說(shuō)明本文在處理大密度比界面問(wèn)題中是合理的。

圖7 低密度流體y-t曲線(xiàn)(算例2)

3.3 非Boussinesq鎖定交換問(wèn)題(算例3)

擋板隔開(kāi)矩形水平管兩側(cè)不同密度的流體(圖8)，當(dāng)瞬間抽去中間擋板(白色區(qū)域)時(shí)，2種流體發(fā)生相對(duì)運(yùn)動(dòng)，其密度比大的情況稱(chēng)之為非Boussinesq鎖定交換問(wèn)題，反之稱(chēng)為Boussinesq鎖定交換問(wèn)題。取矩形管道的長(zhǎng)為1.28 m，寬為0.32 m，擋板左側(cè)為低密度(紅色)流體，右側(cè)為高密度(藍(lán)色)流體(圖8)，各變量見(jiàn)表1。

圖8 初始界面示意圖

表1 擋板左右兩側(cè)參數(shù)表

從圖9可以看出，低密度流體與高密度流體分別沿水平管上、下壁面自右至左和自左至右運(yùn)動(dòng)，形成剪切界面。隨著計(jì)算時(shí)間的增長(zhǎng)，左右相對(duì)運(yùn)動(dòng)的界面存在著切向速度差，出現(xiàn)了Klelvin- Helmholtz界面不穩(wěn)定性。本文通過(guò)使用對(duì)稱(chēng)形式的粒子間壓力，使得界面兩側(cè)壓力符合牛頓第二定律，保證了界面兩側(cè)壓力的連續(xù)性，使得交界面更加光滑，同時(shí)在界面處對(duì)異相粒子施加人工應(yīng)力，避免粒子在異相液體相中的滯留情況。實(shí)驗(yàn)結(jié)果表明(圖9(a)和圖9(b))，本文方法的效果更好，異相流體界面分層更加清晰。

圖9 運(yùn)動(dòng)后界面效果圖

圖10給出輕重流體鋒面與擋板距離/隨時(shí)間變化的曲線(xiàn)。從中可以看出，曲線(xiàn)的變化趨勢(shì)相同，二者近似呈線(xiàn)性變化，與文獻(xiàn)[20]計(jì)算和實(shí)驗(yàn)結(jié)果相符合。其中對(duì)重流體，本文實(shí)驗(yàn)結(jié)果略低于計(jì)算結(jié)果值，初始誤差較小，誤差隨增大而增大并趨于恒定。對(duì)于輕流體，本文實(shí)驗(yàn)結(jié)果與計(jì)算結(jié)果相比，初始時(shí)誤差較大，隨著計(jì)算時(shí)間不斷增大，曲線(xiàn)逐漸趨于統(tǒng)一，誤差趨近于0。說(shuō)明本文方法在計(jì)算大密度比界面問(wèn)題中是合理的。

圖10 輕重流體(x/H)-t曲線(xiàn)

3.4 空氣中圓形液滴形成變化實(shí)驗(yàn)(算例4)

實(shí)驗(yàn)物質(zhì)參數(shù)見(jiàn)表2。為驗(yàn)證多相流大密度比界面處的表面張力計(jì)算公式的有效性，本文進(jìn)行含有表面張力的經(jīng)典實(shí)驗(yàn)——圓形液滴形成變化仿真實(shí)驗(yàn)。初始形態(tài)如圖11所示，其中，黃色粒子表示水粒子，紫色粒子為空氣粒子。

表2 實(shí)驗(yàn)中的物質(zhì)參數(shù)

圖11 空氣中圓形液滴形成變化初始圖

液滴在表面張力的作用下不斷振蕩，由于黏性使得動(dòng)能不斷耗散，最終穩(wěn)定成圓形液滴，如圖12所示。文獻(xiàn)[21]是當(dāng)前模擬該實(shí)驗(yàn)效果最好的方法。將本文SPH計(jì)算結(jié)果與文獻(xiàn)[21]方法進(jìn)行對(duì)比，發(fā)現(xiàn)本文方法產(chǎn)生的液滴更圓，粒子排布更為均勻，在水-氣交界面處更為光滑，從而驗(yàn)證了本文基于密度權(quán)重計(jì)算表面張力公式的有效性。

圖12 圓形液滴形成變化過(guò)程效果圖

本文主要針對(duì)于大密度比兩相流進(jìn)行研究，其中對(duì)稱(chēng)形式的壓力梯度粒子近似公式主要適用于緩慢運(yùn)動(dòng)的兩相流交界面處的壓力梯度計(jì)算，而基于密度權(quán)重的表面張力計(jì)算主要適用于大密度比兩相流在界面處劇烈運(yùn)動(dòng)，從而需要計(jì)算表面張力帶來(lái)的影響的情況。本文尚未考慮是否可以將本文的方法應(yīng)用到三相流甚至更多相流的流體模擬中，同時(shí)基于密度權(quán)重的表面張力計(jì)算主要針對(duì)的是二維情況下的方形液體變化情況，本文的方法是否可以拓展到三維，這些都將是今后研究的重點(diǎn)。

4 結(jié) 論

針對(duì)多相流界面處密度差導(dǎo)致界面壓力與表面張力計(jì)算錯(cuò)誤的問(wèn)題，本文提出了一種新的壓力梯度近似公式和改進(jìn)的界面處人工斥力公式，保證了高密度流體與低密度流體在交界面兩側(cè)的壓力對(duì)稱(chēng)，符合牛頓第三定律，使得多相流交界面更加清晰、光滑、無(wú)穿透。此外，本文還針對(duì)大密度比多相流(例如空氣-水)，給出基于密度權(quán)重的色函數(shù)計(jì)算公式，對(duì)交界面處的表面張力計(jì)算公式進(jìn)行了改進(jìn)；仿真試驗(yàn)結(jié)果表明，本文的表面張力計(jì)算公式在計(jì)算大密度比多相流交界面張力時(shí)，界面更加光滑，模擬效果更好。

[1] 柳有權(quán), 劉學(xué)慧, 朱紅斌, 等. 基于物理的流體模擬動(dòng)畫(huà)綜述[J]. 計(jì)算機(jī)輔助幾何設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2005, 17(12): 2581-2589. LIU Y Q, LIU X H, ZHU H B, et al. Physically based fluid simulation in computer animation[J]. Journal of Computer-Aided Design & Computer Graphics, 2005, 17(12): 2581-2589 (in Chinese).

[2] COLAGROSSI A, LANDRINI M. Numerical simulationof interfacial flows by smoothed particle hydrodynamics[J]. Journal of Computational Physics, 2003, 191(2): 448-475.

[3] MONAGHAN J J, KOCHARYAN A. SPH simulationof multi-phase flow[J]. Computer Physics Communications, 1995, 87(1-2): 225-235.

[4] HU X Y, ADAMS N A. A multi-phase SPH method for macroscopic and mesoscopic flows[J]. Journal of Computational Physics, 2006, 213(2): 844-861.

[5] HU X Y, ADAMS N A. An incompressible multi-phase SPH method[J]. Journal of Computational Physics, 2007, 227(1): 264-278.

[6] SOLENTHALER B, PAJAROLA R. Density contrast SPH interfaces[J]. In Proceedings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, Eurographics Association. New York: IEEE Press, 2008: 211-218.

[7] MONAGHAN J J, RAFIEE A. A simple SPH algorithm for multi-fluid flow with high density ratios[J]. International Journal for Numerical Methods in Fluids, 2013, 71(5): 537-561.

[8] CHEN Z, ZONG Z, LIU M B, et al. An SPH model for multiphase flows with complex interfaces and large density differences[J]. Journal of Computational Physics, 2015, 283: 169-188.

[9] LIND S J, STANSBY P K, ROGERS B D. Incompressible-compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics (SPH)[J]. Journal of Computational Physics, 2015, 309: 129-147.

[10] KRIMI A, REZOUG M, KHELLADI S, et al. Smoothed particle hydrodynamics: a consistent model for interfacial multiphase fluid flow simulations[J]. Journal of Computational Physics, 2018, 358: 53-87.

[11] 楊秋足, 徐緋, 王璐, 等. 一種基于黎曼解處理大密度比多相流SPH的改進(jìn)算法[J]. 力學(xué)學(xué)報(bào), 2019, 51(3): 730-742.YANG Q Z, XU F, WANG L, et al. An improved SPH algorithm for large density ratios multiphase flows based on riemann solution[J]. Chinese Journal of Thoretical and Applied Mechanics, 2019, 51(3): 730-742 (in Chinese).

[12] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astrophysical Journal, 1977, 8(12): 1013-1024.

[13] GINGOLD R A, MONAGHAN J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society 1977, 181(3): 375-389.

[14] LIBERSKY L D, PETSCHEK A G. Smooth particle hydrodynamics with strength of materials[EB/OL]. [2020-03-10]. https://link.springer.com/chapter/10.1007% 2F3-540-54960-9_58.

[15] MONAGHAN J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics, 1994, 110(2): 399-406.

[16] EDMONF Y M, SHAO S D. Simulation of near-shore solitary wave mechanics by an incompressible SPH method[J]. Applied Ocean Research, 2002, 24(5): 275-286.

[17] MONAGHAN J J, GINGOLD R A. Shock simulation by the particle method SPH[J]. Journal of Computational Physics, 1983, 52(2): 374-389.

[18] ADAMI S, HU X Y, ADAMS N A. A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation[J]. Journal of Computational Physics, 2010, 229(13): 5011-5021.

[19] LAWRIE A G W, DALZIEL S B. Turbulent diffusion in tall tubes. I. Models for Rayleigh-Taylor instability[J]. Physics of Fluids, 2011, 23(8): 3097.

[20] LOWE R J, ROTTMAN J W, LINDEN P F. The non-Boussinesq lock-exchange problem, Part 1 theory and experiments[J]. Journal of Fluid Mechanics, 2005, 537: 101-124.

[21] 沈雁鳴, 陳堅(jiān)強(qiáng). SPH方法對(duì)氣液兩相流自由界面運(yùn)動(dòng)的追蹤模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2012(2): 23-27, 34. SHEN Y M, CHEN J Q. Tracking simulation of free interface motion of gas-liquid two-phase flow by SPH method[J]. Acta Aerodynamica Sinica, 2012(2): 23-27, 34 (in Chinese).

Simulation of free interface motion of two-phase flow based on modified pressure and surface tension calculation

ZHU Xiao-lin, ZHOU Yun-ruo, HE Hong-hong

(School of Mathematics, Hefei University of Technology, Hefei Anhui230009, China)

In the simulation of multiphase flow using the smoothed particle hydrodynamics (SPH) method, the density of the multiphase flow at the interface was discontinuous, the pressure calculation of the particle interface was erroneous, and problems occurred, such as pressure oscillation at the interface and interface fracture. In response to the above problems, a new pressure gradient approximation formula and an improved artificial repulsion formula at the interface were proposed to make the interface clearer and smoother without penetration, thus producing better simulation results. In addition, by giving a color function calculation formula based on density weights, the formula was improved for the calculation of surface tension at the interface of multiphase flow with high density ratio, leading to the smoother density transition of multiphase flow interface and better simulation effect. Finally, through three simulation experiments, such as dam break, Rayleigh-Taylor interface instability, and non-Boussinesq locked exchange problem, the particle distribution map of the interface and the distance of the interface front at different times were obtained. The results can verify the rationality of the new pressure gradient approximation formula and the artificial repulsion formula at the interface. Through the simulation experiment of droplet formation in the air, the particle change diagram formed by the circular droplet was obtained, and the result can corroborate the effectiveness of the improved surface tension calculation method in this paper.

smoothed particle hydrodynamics method; multi-phase flow; pressure gradient particle approximation; artificial repulsion; the color function; surface Tension

TP 391

10.11996/JG.j.2095-302X.2020060970

A

2095-302X(2020)06-0970-10

2020-06-22；

2020-07-28

22 June，2020；

28 July，2020

國(guó)家自然科學(xué)基金項(xiàng)目(61272024)

National Natural Science Foundation of China (61272024)

朱曉臨(1964-)，男，安徽池州人，教授，博士，碩士生導(dǎo)師。主要研究方向?yàn)閿?shù)值逼近、計(jì)算機(jī)圖形學(xué)、CAGD、圖形圖像處理等。 E-mail：zxl_hfut@126.com

ZHU Xiao-lin (1964-), male, professor, Ph.D. His main research interests cover numerical approximation, computer graphics, CAGD and graphic image processing. E-mail：zxl_hfut@126.com