劉啟新,孫中國(guó),段壯,席光

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,710049,西安)

無(wú)網(wǎng)格粒子法,如移動(dòng)粒子半隱式(MPS)法[1]和光滑粒子動(dòng)力學(xué)(SPH)法[2]是一類(lèi)建立在拉格朗日框架下的計(jì)算流體動(dòng)力學(xué)數(shù)值方法。相較于歐拉框架下的網(wǎng)格法,無(wú)網(wǎng)格粒子法在具有高度非線性、充分發(fā)展自由表面和大變形等特征的流動(dòng)問(wèn)題的研究上具有更大的優(yōu)勢(shì)。因此,這類(lèi)方法在海洋工程[3-4]、化工過(guò)程[5-6]、核工程[7-8],多相流[9-10]等領(lǐng)域都已得到非常廣泛應(yīng)用。隨著物理模型、數(shù)學(xué)算法的不斷完善以及計(jì)算機(jī)計(jì)算能力的快速發(fā)展,無(wú)網(wǎng)格粒子法的應(yīng)用正在向更廣泛的領(lǐng)域擴(kuò)展。

在前處理環(huán)節(jié)中,粒子法通常使用間距相等的計(jì)算節(jié)點(diǎn)離散空間中的連續(xù)計(jì)算域,這些計(jì)算節(jié)點(diǎn)被形象地表述為有統(tǒng)一固定體積的光滑粒子。離散代數(shù)方程對(duì)每一個(gè)粒子所攜帶的物理量信息(空間坐標(biāo)、速度、壓力、溫度等)進(jìn)行求解,該過(guò)程不需要粒子間具有明確的拓?fù)浣Y(jié)構(gòu),只需要確定粒子間實(shí)時(shí)的位置關(guān)系。因此,相較于網(wǎng)格法,粒子法的前處理環(huán)節(jié)一直以來(lái)更為簡(jiǎn)單,無(wú)論是二維還是三維問(wèn)題,都可以采用矩形陣列的初始粒子分布。

但是,隨著粒子法在工業(yè)領(lǐng)域的日益發(fā)展,所要離散的流體域構(gòu)型越來(lái)越復(fù)雜,在應(yīng)對(duì)彎曲或傾斜壁面時(shí),矩形陣列分布的粒子難以保證貼體均勻的分布,而不均勻分布的粒子往往更容易引起數(shù)值振蕩[11]。因缺乏有效的前處理粒子布置方案,若想要得到貼體、均勻、各向同性的粒子分布,就需要基于拉格朗日粒子法的特性,按照進(jìn)口邊界條件使粒子流進(jìn)系統(tǒng),逐漸充滿復(fù)雜流體域[12-14]。這就相當(dāng)于在計(jì)算時(shí)要經(jīng)過(guò)完整的非定常啟動(dòng)過(guò)程,而這個(gè)過(guò)程對(duì)于數(shù)值研究通常不是必須的,這就增加了計(jì)算需要的資源消耗和研究周期。因此,需要為粒子法工業(yè)研究提出更規(guī)范的前處理方法。

目前,粒子布置法的研究主要關(guān)注粒子對(duì)固體邊界構(gòu)型的取樣。李帝辰等[15]提出了將二維問(wèn)題降為一維問(wèn)題的粒子構(gòu)型法,該方法對(duì)復(fù)雜邊界進(jìn)行均勻的表面粒子分布之后再進(jìn)行內(nèi)部填充。Domínguez等[16]開(kāi)發(fā)了一款名為GenCase的粒子法前處理軟件,該軟件通過(guò)網(wǎng)格三維模型的節(jié)點(diǎn)建立粒子的初場(chǎng)布置,能夠?qū)崿F(xiàn)各種復(fù)雜三維構(gòu)型的粒子初場(chǎng)布置,但該軟件不能保證邊界粒子對(duì)光滑曲面的貼體性。Akinci等[17]考慮了固體邊界對(duì)流體粒子物理量的計(jì)算,使用單層粒子采樣復(fù)雜固壁;Ihmsen等[18]針對(duì)粒子法的流固耦合問(wèn)題提出了一種預(yù)測(cè)-修正的固體邊界處理方法,在考慮固體邊界的壓力情況下計(jì)算,但Akinci和Ihmsen的研究?jī)?nèi)容并不涉及初始布置。倪國(guó)喜等[19]和Diehl等[20]基于Voronoi方法的網(wǎng)格重構(gòu)法分別建立了二維和三維的均勻各向同性粒子初場(chǎng)布置算法;Zhu等[21]基于level-set方法構(gòu)建了一種復(fù)雜固體壁面粒子離散的方案,但以上方法都需要引入全新的數(shù)值方法,一定程度上增加了粒子法求解器前處理模塊編寫(xiě)的復(fù)雜度和難度。

考慮到目前算法耦合在粒子法工業(yè)應(yīng)用中越發(fā)廣泛,越來(lái)越多的工程問(wèn)題開(kāi)始采用粒子法和網(wǎng)格法相耦合的方式來(lái)研究復(fù)雜的流固耦合問(wèn)題,如流體機(jī)械、過(guò)程裝備、海洋船舶、核反應(yīng)設(shè)備等問(wèn)題,其中粒子法模擬流體部分,網(wǎng)格法模擬固體部分。因這些問(wèn)題的流體計(jì)算域具有復(fù)雜的幾何結(jié)構(gòu),因而亟需更合理、更高效的粒子初場(chǎng)布置方案。

基于以上,本文嘗試借助網(wǎng)格邊界,將基于該網(wǎng)格邊界的數(shù)值積算法引入基于濃度梯度的粒子遷移模型中,建立起一套簡(jiǎn)潔、靈活的粒子初場(chǎng)布置的優(yōu)化方法,使矩形陣列分布的初場(chǎng)粒子進(jìn)行空間位置的自適應(yīng)調(diào)整,以提高粒子分布的均勻各向同性和固體邊界附近的流體粒子分布對(duì)該邊界的貼體性。此外,粒子法中的自由表面粒子表征液體和外界大氣環(huán)境的交界面,在壓力計(jì)算中扮演著自由表面壓力邊界條件的重要角色?；谝话懔Ｗ臃ńy(tǒng)一計(jì)算節(jié)點(diǎn)尺度的特點(diǎn),原本僅有分子自由程尺度的氣液交界面在粒子法中具有和粒子直徑同等尺度的“厚度”。因此,要求計(jì)算域表面粒子分布光滑、連續(xù),符合設(shè)定的自由邊界,才能保證自由邊界附近的壓力計(jì)算的正確性。例如,液滴的初始布置的理想邊界形狀為圓或球體,駐波的初始邊界的波形為三角函數(shù)曲線等。本文展示了所提算法在自由邊界的初場(chǎng)布置上的應(yīng)用,將初始固體表面(固-液交界面)和初始自由表面(氣-液交界面)均定義為理論界面,所提方法可以有效提高初始布置中界面粒子分布的精度。

1 數(shù)值算法

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

本文以移動(dòng)粒子半隱式(MPS)法為例,介紹本文提出的初始布置優(yōu)化算法在粒子法中的應(yīng)用。MPS法的控制方程為不可壓縮NS方程,表示為

(1)

(2)

式中:u為速度矢量;t為時(shí)間;ρ為密度;P為壓力;ν為運(yùn)動(dòng)黏度;g為重力加速度;fs為表面張力。

粒子間的相互作用模型是基于核函數(shù)的加權(quán)平均思想建立起來(lái)的,核函數(shù)為

(3)

re=2.5l0

(4)

式中:|rij|為i粒子和j粒子間的位置矢量的模;re為影響域半徑;l0為粒子間距離(或每一個(gè)粒子的直徑)。

求解NS方程時(shí)的微分算子,包括梯度算子、散度算子和Laplace算子,如下式所示

(5)

(6)

(7)

式中:d表示維度;φ代表壓力、溫度等物理量;n0表示初始粒子數(shù)密度,其定義如下

(8)

關(guān)于MPS法的更細(xì)節(jié)的算法,可參考文獻(xiàn)[1]。

1.2 粒子遷移模型

基于濃度梯度的粒子遷移模型最早由Lind等[22]提出,其目標(biāo)是改善高精度粒子法中流體粒子因嚴(yán)格跟隨跡線運(yùn)動(dòng)所產(chǎn)生的粒子分布不均的問(wèn)題。無(wú)論是對(duì)于全顯式求解的SPH法,還是半隱式求解的MPS法,粒子分布不均的問(wèn)題都會(huì)破壞流場(chǎng)求解的穩(wěn)定性,導(dǎo)致誤差急速累積,致使計(jì)算在短時(shí)間內(nèi)崩潰。因此,在粒子法向提高精度的目標(biāo)不斷發(fā)展的過(guò)程中,粒子遷移模型已逐漸成為求解過(guò)程中不可或缺的模型。

到目前為止,最廣泛采用的粒子遷移模型依然遵循濃度梯度驅(qū)動(dòng)粒子移動(dòng)的思想,如式(9)所示[23]

(9)

式中:δri表示i粒子的遷移矢量;wij為核函數(shù)。

當(dāng)流場(chǎng)中的粒子分布不均勻時(shí),局部粒子周?chē)鷷?huì)產(chǎn)生粒子數(shù)密度梯度,上述粒子遷移模型便可以驅(qū)使粒子從高粒子數(shù)密度區(qū)域向低粒子數(shù)密度區(qū)域移動(dòng),從而使粒子的分布回歸均勻狀態(tài)。本文將這一思想應(yīng)用在初始布置中,通過(guò)粒子遷移模型作用下粒子會(huì)自發(fā)形成均勻的空間分布這一特點(diǎn),對(duì)矩形陣列的初場(chǎng)粒子進(jìn)行空間位置上的重構(gòu),以改善邊界粒子分布的不連續(xù)性,提高對(duì)光滑連續(xù)邊界的貼體性。

1.3 基于邊界網(wǎng)格的數(shù)值積分粒子遷移模型

Zheng等[24]最早將基于面積微元的面積分算法引入到了EMPS法中,代替了原本的固體邊界粒子,實(shí)現(xiàn)了EMPS法中的無(wú)穿透固體壁面邊界條件。本文將面積分算法引入到粒子遷移模型中,建立理論邊界及外場(chǎng)一定范圍內(nèi)的網(wǎng)格模型,即邊界網(wǎng)格,如圖1所示。通過(guò)對(duì)邊界網(wǎng)格建立面積分算法,實(shí)現(xiàn)粒子的均勻分布。

圖1 流體粒子-邊界網(wǎng)格示意圖Fig.1 Schematic diagram of fluid particle-wall meshes

基于邊界網(wǎng)格的粒子數(shù)密度計(jì)算公式如下

(10)

式中:ΔAk表示網(wǎng)格微元k的面積,即圖1中網(wǎng)格部分的每一個(gè)網(wǎng)格微元的面積;xk表示網(wǎng)格微元k的形心坐標(biāo),即圖1中紅色的節(jié)點(diǎn)坐標(biāo);f為流體粒子;m為圖1中的網(wǎng)格。

目前廣泛用于劃分網(wǎng)格的商業(yè)軟件,如ICEM、Gambit等,均可提供網(wǎng)格節(jié)點(diǎn)的坐標(biāo)信息。需要說(shuō)明兩點(diǎn):第一,該算法對(duì)網(wǎng)格法模擬中通常會(huì)關(guān)注的正交性等網(wǎng)格質(zhì)量評(píng)價(jià)沒(méi)有要求;第二,結(jié)構(gòu)化矩形網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格都可以用于該算法。

圖2所示為3種邊界網(wǎng)格。本文通過(guò)計(jì)算流體域的粒子數(shù)密度,驗(yàn)證邊界網(wǎng)格的粒子數(shù)密度計(jì)算公式的正確性。3種邊界網(wǎng)格分別為:①正方形粗邊界網(wǎng)格,網(wǎng)格單元的邊長(zhǎng)與粒子間間距相等;②三角形邊界網(wǎng)格,網(wǎng)格單元的直角邊長(zhǎng)與粒子間間距相等;③正方形細(xì)邊界網(wǎng)格,網(wǎng)格單元的邊長(zhǎng)為粒子間間距的1/2。計(jì)算結(jié)果顯示,在使用3種邊界網(wǎng)格時(shí),計(jì)算得到的粒子數(shù)密度在整個(gè)計(jì)算域中的分布都是一致的,即內(nèi)部粒子因其支撐域的完整性可以達(dá)到初始粒子數(shù)密度n0,n0的值是核函數(shù)和維度的函數(shù),本文中的n0約為2.48;而自由邊界粒子支撐域內(nèi)沒(méi)有足夠粒子,因此粒子數(shù)密度小于n0。粒子數(shù)密度的計(jì)算證明了基于邊界網(wǎng)格的核函數(shù)積分的正確性。

(a)正方形粗邊界網(wǎng)格

(b)三角形邊界網(wǎng)格

(c)正方形細(xì)邊界網(wǎng)格

1.4 基于邊界網(wǎng)格的粒子遷移模型

基于以上結(jié)果,可以得到在使用邊界網(wǎng)格時(shí)的粒子遷移模型,如式(11)所示。

(11)

式中:系數(shù)A是與維度d、粒子直徑l0、粒子數(shù)密度n0有關(guān)的量,其計(jì)算公式如下

(12)

i粒子的遷移量分為兩部分,第一部分是式(11)等號(hào)右端的第一項(xiàng),代表粒子i影響域內(nèi)的其他粒子對(duì)該粒子的遷移量所做出的貢獻(xiàn);第二部分是等號(hào)右端的第二項(xiàng),代表邊界網(wǎng)格對(duì)粒子i的遷移量做出的貢獻(xiàn)。兩部分貢獻(xiàn)線性相加之后即為當(dāng)前i粒子的空間遷移矢量值。

式(12)中的α值用以控制遷移矢量的模的大小,模量越大,單步內(nèi)粒子移動(dòng)的位置就越長(zhǎng)。合適的α值可以同時(shí)保證粒子分布重構(gòu)過(guò)程的高效和穩(wěn)定性。

經(jīng)過(guò)若干步的粒子遷移計(jì)算之后,流場(chǎng)中的粒子會(huì)自發(fā)向延粒子數(shù)密度梯度的反方向運(yùn)動(dòng),最終發(fā)展為均勻各向同性的分布。這種分布也是不可壓縮粒子法計(jì)算中節(jié)點(diǎn)的理想分布。

1.5 CFL條件

本文所提出的算法不包含時(shí)間項(xiàng),粒子重構(gòu)過(guò)程中每一步的移動(dòng)尺度由式(12)中的α控制。綜合考慮粒子遷移過(guò)程能夠穩(wěn)定且高效地推進(jìn),需要?jiǎng)討B(tài)地調(diào)整遷移模型中的系數(shù)α,其設(shè)定依據(jù)以下條件

|δri|max≤0.1l0

(13)

2 結(jié)果與討論

2.1 二維粒子初場(chǎng)布置的重構(gòu)及優(yōu)化過(guò)程

以橢圓計(jì)算域?yàn)槔?進(jìn)行初始粒子分布的各向同性重構(gòu)。橢圓形邊界具有明確的函數(shù)表達(dá),但難以使用圓形粒子布置的方法在極坐標(biāo)下進(jìn)行邊界粒子逐個(gè)配點(diǎn),此外,橢圓內(nèi)部的粒子也難以進(jìn)行均勻各向同性的逐個(gè)配點(diǎn)。

在布置均勻各向同性的橢圓形粒子分布時(shí),以預(yù)先建立的橢圓形網(wǎng)格為輔助計(jì)算節(jié)點(diǎn),如圖3所示。橢圓結(jié)構(gòu)的長(zhǎng)軸長(zhǎng)度與短軸長(zhǎng)度分別為2 m和1 m,主網(wǎng)格尺寸為0.02 m。網(wǎng)格域的寬度采取粒子的影響域半徑re,其數(shù)據(jù)需要包含每一個(gè)網(wǎng)格微元k的形心xk和面積ΔAk。

圖3 橢圓形網(wǎng)格信息Fig.3 The information of the elliptic meshes

之后,建立矩形陣列分布的流場(chǎng)節(jié)點(diǎn),即流體粒子,如圖4(a)所示。粒子直徑為0.1 m,粒子數(shù)為586。矩形陣列的粒子離散域的階梯狀非連續(xù)邊界偏離了橢圓形邊界。因此,需基于方程(11)對(duì)粒子的空間坐標(biāo)進(jìn)行調(diào)整,減小邊界的離散誤差,結(jié)果如圖4(b)所示。

(a)矩形陣列粒子分布

(b)粒子重構(gòu)后的各向同性均勻分布

distribution within the elliptic region

對(duì)重構(gòu)過(guò)程前后兩種粒子分布的離散誤差進(jìn)行量化對(duì)比,結(jié)果如圖5所示。對(duì)比橢圓形計(jì)算域中的兩種粒子分布的粒子數(shù)密度如圖6所示,可以看出,重構(gòu)后的粒子數(shù)密度更接近理論值。從以上兩項(xiàng)定量對(duì)比可以看出,相比起陣列布置中的表面粒子,重構(gòu)后的表面粒子與理論邊界的偏差值和粒子數(shù)密度差整體降低,重構(gòu)后的最大偏差較重構(gòu)前大約減小了81%,粒子數(shù)密度減小了約60%,這表明重構(gòu)后的表面粒子分布更加貼合理論邊界。

圖5 橢圓形計(jì)算域重構(gòu)過(guò)程前后表面粒子與理論邊界的偏差Fig.5 The deviations between the surface particles and the theoretical boundary before and after the reconstruction in the elliptic domain

圖6 橢圓計(jì)算域重構(gòu)過(guò)程前后表面粒子的粒子數(shù)密度Fig.6 Particle number density of surface particles before and after the reconstruction in elliptic domain

將以上兩種橢圓分布視為初始布置,將橢圓形邊界視為自由表面,使用移動(dòng)粒子半隱式法[25],考察在保守力場(chǎng)中初始時(shí)刻的壓力

P=-Ωr2

(14)

計(jì)算域初始速度定義為

(15)

式中:Ω和δ(0)為定義保守力場(chǎng)大小的兩個(gè)常數(shù),值分別為0.4和1.2。圖7對(duì)兩種粒子分布下的物理量場(chǎng)計(jì)算進(jìn)行了對(duì)比,并對(duì)粒子所離散的物理量場(chǎng)進(jìn)行連續(xù)化處理,直觀地展示了理想的粒子分布對(duì)物理量場(chǎng)求解帶來(lái)的改善效果。圖7(a)中,呈階梯狀的階躍邊界粒子(壓力為零粒子)所構(gòu)成的自由邊界會(huì)造成內(nèi)部壓力場(chǎng)求解失真,而圖7(b)中的連續(xù)光滑邊界則提供了更合理的自由表面壓力邊界條件,基本避免了粒子法自由邊界“厚度”對(duì)計(jì)算造成的影響。

(a)矩形陣列粒子分布

(b)粒子重構(gòu)后的各向同性均勻分布

2.2 二維轉(zhuǎn)子泵內(nèi)流場(chǎng)的粒子初始布置

二維轉(zhuǎn)子泵轉(zhuǎn)子與轉(zhuǎn)子之間、轉(zhuǎn)子和泵腔間的嚙合使泵腔內(nèi)空間結(jié)構(gòu)復(fù)雜,尤其在嚙合點(diǎn)附近的狹縫區(qū)域,粒子難以用配點(diǎn)法進(jìn)行手動(dòng)布置。

圖8展示了轉(zhuǎn)子泵固體部件的網(wǎng)格結(jié)構(gòu),網(wǎng)格形態(tài)均為三角形,網(wǎng)格主尺寸約為5 m。粒子的初場(chǎng)布置如圖9(a)所示,其中,粒子直徑為2.5 m,粒子數(shù)為7 574。

圖8 轉(zhuǎn)子泵網(wǎng)格Fig.8 The rotor pump grid

(a)矩形陣列粒子分布

(b)粒子重構(gòu)后的各向同性均勻分布

由圖9(a)的放大圖可以看到,矩形陣列的粒子難以均勻地填充其嚙合點(diǎn)附近的狹縫,同時(shí),固壁邊界附近的粒子因粒子在空間位置中的截?cái)鄬?dǎo)致當(dāng)?shù)亓Ｗ臃植寂c網(wǎng)格的拓?fù)浣Y(jié)構(gòu)缺乏連續(xù)性,因此會(huì)產(chǎn)生較嚴(yán)重的粒子數(shù)密度噪點(diǎn)。圖9(b)是經(jīng)過(guò)了粒子分布重構(gòu)之后的結(jié)果。從局部放大圖中可以看到,經(jīng)過(guò)自適應(yīng)調(diào)整后的粒子分布保證了均勻各向同性,與固體邊界之間具有更好的貼體性和空間連續(xù)性,并且更符合嚙合點(diǎn)附近不斷收縮的區(qū)域形狀。

圖10同樣對(duì)比了重構(gòu)前后表面粒子的粒子數(shù)密度。可以看出,重構(gòu)后的固壁邊界附近的粒子數(shù)密度在坐標(biāo)系中的分布更加集中,且更加接近理論粒子數(shù)密度,數(shù)據(jù)波動(dòng)更小,方差更低。

圖10 轉(zhuǎn)子泵重構(gòu)前后表面粒子的粒子數(shù)密度Fig.10 Particle number density of surface particles before and after the reconstruction of particle distribution in the rotor pump

2.3 三維復(fù)雜空間中的粒子初場(chǎng)布置

本節(jié)繼續(xù)考察本研究所提出的粒子初始布置重構(gòu)算法在復(fù)雜三維計(jì)算域中的應(yīng)用。如圖11所示,在由立方體外邊界和五角星內(nèi)邊界組合而成的復(fù)雜結(jié)構(gòu)空腔中,需要流體粒子均勻充滿該空腔。為便于展示,圖中以三維實(shí)心粒子表示流體粒子與網(wǎng)格微元的形心坐標(biāo),并隱藏網(wǎng)格在Z方向上表面的形心坐標(biāo)。同樣使用網(wǎng)格將固體部分離散,流體粒子被均勻地陣列在空腔中。算例中,粒子直徑為0.01 m,網(wǎng)格主尺寸為0.007 m,粒子數(shù)為28 784。

圖11 復(fù)雜形狀腔體的初始粒子布置和網(wǎng)格形心Fig.11 Initial particle distribution and surrounding mesh center of a complex cavity

均勻化前后的粒子分布如圖12所示,可以看出,矩陣陣列的粒子在五角星的斜邊上形成的階梯狀不連續(xù)分布在重構(gòu)后得以改善,更加貼合傾斜直邊。同時(shí),在五角星的尖角處,粒子環(huán)繞尖角的分布也更加合理。這種銳角結(jié)構(gòu)中的粒子分布在粒子法的布置中屬于較為困難的例子,而本文提出的重構(gòu)方案,使這種特殊結(jié)構(gòu)的粒子布置得到了很好的處理。此外,從物理場(chǎng)云圖可以看出,整場(chǎng)的粒子數(shù)密度在重構(gòu)后也更加均勻,這有助于有效抑制計(jì)算初始時(shí)刻的壓力振蕩,保證了計(jì)算的穩(wěn)定性。

量化對(duì)比了重構(gòu)前和重構(gòu)后邊界附近的粒子數(shù)密度,結(jié)果如圖13所示。重構(gòu)后的粒子數(shù)密度更集中地分布于理論粒子數(shù)密度的區(qū)域,最大誤差相較于重構(gòu)前減小了約51%,反映了更均勻的粒子數(shù)密度,能夠有效提高粒子法的水動(dòng)力學(xué)模擬質(zhì)量。

2.4 三維球體液滴的粒子初始布置及浸潤(rùn)計(jì)算

2.4.1 三維液滴的初始布置

三維球體是液滴動(dòng)力學(xué)問(wèn)題的研究中廣泛存在的基本結(jié)構(gòu),然而其離散粒子的初場(chǎng)布置很難還原球體表面。就作者所知,目前廣泛采用的三維液滴初場(chǎng)布置方案依然采用截取球心的半徑范圍內(nèi)均勻陣列布置的流體粒子。當(dāng)研究復(fù)雜的液滴動(dòng)力學(xué)問(wèn)題時(shí),如單液滴浸潤(rùn)固體壁面、多液滴間的聚并和碰撞等,在這種不理想的分布下,階躍、不連續(xù)的自由表面會(huì)造成表面曲率和法向量的嚴(yán)重誤差,進(jìn)而影響表面張力和壓力的正確計(jì)算,導(dǎo)致表面近表面粒子的壓力振蕩,引起計(jì)算誤差,或在液滴表面產(chǎn)生非物理的流動(dòng);此外,液滴內(nèi)部的各向異性分布也會(huì)造成非物理的壓力噪點(diǎn),致使液滴形態(tài)的模擬受到影響。因此,光滑、連續(xù)的二維圓形或三維球體表面粒子分布對(duì)提高液滴動(dòng)力學(xué)研究的嚴(yán)謹(jǐn)性和規(guī)范性具有重要意義。

目前,在絕大多數(shù)研究三維液滴動(dòng)力學(xué)特性的粒子法模擬中,液滴的初始布置都采用圖14(a)的布置方案。在這種初始分布下的液滴形態(tài)的模擬會(huì)因液滴內(nèi)部分布在同一平面內(nèi)的粒子沿平面方向的擠壓而失真,而經(jīng)過(guò)優(yōu)化的各向同性粒子分布在計(jì)算過(guò)程中能夠保持更合理的形態(tài),如圖14(b)所示。

(a)三維球體的矩形陣列粒分布

(b)三維球體的均勻各向同性粒子分布

2.4.2 三維液滴的浸潤(rùn)計(jì)算

圖15展示了兩種粒子分布下5、25 ms的浸潤(rùn)過(guò)程。5 ms時(shí),液滴在固液界面張力的作用下開(kāi)始浸潤(rùn)固壁。未重構(gòu)的液滴在計(jì)算初期,其表面粒子因沿自由表面法向的表面張力和壓力梯度之間的相互作用而發(fā)生非物理的振蕩,因此其自由表面相比起重構(gòu)后的初始布置更加粗糙。此外,未重構(gòu)的粒子內(nèi)因粒子各向異性分布,內(nèi)部同一平面內(nèi)的粒子沿平面方向的擠壓,破壞液滴模擬過(guò)程中的液滴圓度。在25 ms時(shí),原始均勻陣列的液體形態(tài)已明顯失真,呈現(xiàn)出近四邊形,而優(yōu)化后的各向同性初始分布的液滴形態(tài)更符合實(shí)際。

(a)5 ms

(a)25 ms

量化的圓度對(duì)比見(jiàn)圖16。其中,圓度由XY平面內(nèi)液滴最大截面的外切圓直徑lo與內(nèi)切圓直徑li之比表示。圖中,以原始矩形陣列分布的粒子作為初場(chǎng)布置,液滴在浸潤(rùn)過(guò)程中,XY平面內(nèi)的最大截面的圓度在理論圓度(lo/li=1)附近存在較大波動(dòng),而基于重構(gòu)后的初場(chǎng)布置計(jì)算得到的結(jié)果更接近理論圓度,最大偏差較未重構(gòu)時(shí)減小了約78.6%。該形態(tài)學(xué)上的改善對(duì)于提高小韋伯?dāng)?shù)液滴動(dòng)力學(xué)問(wèn)題的粒子法模擬精度具有重要意義。

圖16 傳統(tǒng)陣列初始布置和優(yōu)化初始布置下浸潤(rùn)過(guò)程的圓度對(duì)比Fig.16 The comparison of roundness during the wetting process between the original initial particle distribution and the reconstituted initial particle distribution

3 結(jié) 論

考慮到粒子法在復(fù)雜工程問(wèn)題中越來(lái)越廣泛的應(yīng)用,為了優(yōu)化具有復(fù)雜邊界的粒子初場(chǎng)布置,本文將粒子法中廣泛采用的粒子遷移模型與粒子法離散模型的數(shù)值積分形式相結(jié)合,建立起一種面向粒子法求解器前處理過(guò)程中粒子分布重構(gòu)的方法。該方法通過(guò)粒子遷移模型,使矩形陣列排布的粒子通過(guò)自適應(yīng)運(yùn)動(dòng),距離設(shè)定邊界(固壁或自由表面)最近的粒子形成更貼合邊界的粒子分布,而邊界附近的粒子也能保證密度均勻、各向同性的分布。

通過(guò)二維橢圓形計(jì)算域、轉(zhuǎn)子泵泵腔、三維復(fù)雜腔體、三維球體等初始布置,驗(yàn)證了算法在含有曲線面、狹縫、銳角等形狀的邊界中具有通用性和有效性。重構(gòu)后的表面粒子所表示的計(jì)算域邊界與理論邊界的偏差值比未重構(gòu)時(shí)減小約81%,且重構(gòu)后的表面粒子數(shù)的噪點(diǎn)更低,整體更接近理論值,其中,二維情況下誤差可減小約60%,三維情況下誤差可減小51%。在液滴浸潤(rùn)算例中,重構(gòu)后的三維液滴避免了各向異性粒子分布導(dǎo)致的形態(tài)失真問(wèn)題,因此在浸潤(rùn)過(guò)程中圓度與理論值符合程度較高,比未重構(gòu)時(shí)誤差減小約78.6%,得到了準(zhǔn)確的液滴形態(tài)。