韓亞偉 強(qiáng)洪夫 趙玖玲 高巍然

(第二炮兵工程大學(xué)601室,西安 710025)

(2012年8月10日收到;2012年9月16日收到修改稿)

1 引言

近年來,光滑粒子流體動(dòng)力學(xué)(SPH)方法作為一種純Lagrangian型的無網(wǎng)格算法得到了很快發(fā)展,該算法最早是由Lucy[1]與Gingold[2]分別提出,最初用于解決天體物理中無邊界三維流動(dòng)問題,現(xiàn)已推廣應(yīng)用到計(jì)算力學(xué)諸多領(lǐng)域的數(shù)值仿真之中.該算法的特點(diǎn)是在計(jì)算中不需要生成任何輔助網(wǎng)格,因此在處理大變形問題時(shí)不存在Lagrangian網(wǎng)格算法中網(wǎng)格纏繞的技術(shù)瓶頸,理論上可以處理任意的變形問題;并且由于算法的Lagrangian特性,SPH離散粒子自然地追蹤流體物質(zhì)的運(yùn)動(dòng)及邊界的變化,不需要像Eluerian網(wǎng)格算法中采取復(fù)雜的界面追蹤技術(shù),不存在Eluerian網(wǎng)格算法中因數(shù)值耗散造成界面追蹤精度的降低.由于SPH算法的以上特點(diǎn),使得它特別適合模擬自由表面、運(yùn)動(dòng)界面及大變形問題[3?6].

SPH方法發(fā)展還尚未成熟,尤其是固壁邊界條件的施加方法,一直是困擾SPH方法發(fā)展的難點(diǎn)之一.早期Monaghan[7]模擬高速自由液面流動(dòng)時(shí)引入了Lennard-Jones對(duì)勢(shì)的強(qiáng)排斥力,其中固壁邊界用虛粒子來描述,該模型中虛粒子對(duì)鄰近邊界的粒子強(qiáng)制施加排斥力,從而阻止鄰近邊界的粒子非物理穿透邊界.但是該方法動(dòng)量不守恒,對(duì)鄰近邊界的粒子的描述與實(shí)際不符.Libersky等[8]和Randles等[9]首次提出鏡像粒子的方法來施加固壁邊界,在邊界之外布置關(guān)于邊界對(duì)稱的鏡像粒子,該粒子與真實(shí)粒子密度相同,速度相反,并參與計(jì)算,該方法守恒性好但鏡像粒子生成復(fù)雜,無法處理復(fù)雜壁面邊界的問題.Gotoh等[10]建議固壁邊界采用邊界粒子平衡內(nèi)部壓強(qiáng)來防止粒子穿透邊界,界面粒子采用多層鏡像粒子來實(shí)現(xiàn)滑移邊界條件,但仍難以處理復(fù)雜邊界問題.文獻(xiàn)[11,12]采用兩種類型的虛粒子來處理固定邊界條件,即一種虛粒子設(shè)置在固定邊界上,與Monaghan[7]所用的相似,另一種虛粒子分布在邊界的周圍,這與Libersky等[8]的鏡像粒子相似,即同時(shí)使用兩種虛粒子來阻止實(shí)粒子穿透邊界,但是這一方法也繼承了兩類虛粒子方法的不足.為了能夠處理復(fù)雜壁面邊界問題,近期國(guó)內(nèi)外學(xué)者把研究的重點(diǎn)放在如何更好地施加排斥力上,先后提出了幾種排斥力公式.Rogers等[13]提出了一種排斥力方法,該方法能很好地防止流體穿過壁面,但不能保證動(dòng)量守恒;強(qiáng)洪夫等[14]基于Galerkin形式的SPH方程推導(dǎo)的罰函數(shù)邊界力能較好地保證動(dòng)量守恒,但該方法含有未知參數(shù),需要經(jīng)過試算才能確定其取值,且邊界附近的流體粒子存在運(yùn)動(dòng)不穩(wěn)定問題;Monaghan等[15]和Liu等[16]也分別提出了各自的排斥力施加方法,這兩種方法類似,都能處理復(fù)雜邊界,但都需要將邊界粒子尺寸設(shè)為流體粒子尺寸的一半或更小,從而導(dǎo)致邊界粒子數(shù)目的增加,降低了計(jì)算效率,且這兩種方法不能保證動(dòng)量守恒.

為了既保持排斥力方法能處理復(fù)雜邊界的優(yōu)勢(shì),又能克服上述幾種方法的不足,本文提出一種新的排斥力模型來施加固壁邊界條件.分別通過靜止液柱算例和液柱坍塌算例,驗(yàn)證了各種排斥力方法的動(dòng)量守恒特性及邊界粒子尺寸對(duì)排斥力施加效果的影響;當(dāng)僅采用邊界上單層粒子計(jì)算邊界力時(shí)存在邊界粒子缺失問題,本文采用了Liu等[16]提出的耦合邊界法中的虛粒子處理方法,以提高邊界附近流體粒子的計(jì)算精度,將該法應(yīng)用到容器中液體的靜止算例及潰壩算例.對(duì)容器中液體的靜止算例,計(jì)算得到了容器底部壓強(qiáng)變化過程,通過與Monaghan等[15]與及強(qiáng)洪夫[14]所提出的方法進(jìn)行對(duì)比,本文方法可有效克服邊界附近的流體粒子存在的運(yùn)動(dòng)不穩(wěn)定問題;通過潰壩算例驗(yàn)證了本文方法處理復(fù)雜邊界問題的有效性.

2 SPH的基本方程及離散

2.1 控制方程

拉格朗日描述的流體動(dòng)力學(xué)方程為:

參數(shù)P0可通過P0=100ρ0V/γ獲得,其中ρ0為流體初始密度,γ=7,Vmax為流體的最大速度.為保證流體的弱可壓縮性,一般取聲速c=10Vmax.

2.2 流體動(dòng)力學(xué)方程的離散

在SPH方法中,連續(xù)的流場(chǎng)離散成為一系列相互作用的粒子,通過核函數(shù)估計(jì)技術(shù)在這些粒子上離散控制方程組,得到一組描述各粒子物理量隨時(shí)間變化的常微分方程組,即SPH基本方程組,再對(duì)這組方程采用相應(yīng)的常微分方程組求解方法來推進(jìn)時(shí)間進(jìn)程的求解.

與耗散粒子動(dòng)力學(xué)方法類似[17],求解域?被離散為一系列質(zhì)點(diǎn),其中每一個(gè)質(zhì)點(diǎn)i記錄了流場(chǎng)在該點(diǎn)處的物理量,如密度ρi,速度vi,位置xi等,并且給予一個(gè)統(tǒng)計(jì)意義的質(zhì)量mi,每一個(gè)粒子同時(shí)擁有一個(gè)影響域S,影響域內(nèi)所包含的其他SPH粒子 j為該粒子的鄰近粒子,對(duì)于系統(tǒng)中任意點(diǎn)x處的物理量及其空間導(dǎo)數(shù)通過核函數(shù)估計(jì)得到,詳見文獻(xiàn)[3,4].

這樣,流體運(yùn)動(dòng)的控制方程(1)—(3)可離散為

2.3 時(shí)間積分

形成虛粒子,從而導(dǎo)致其缺乏通用性,很難推廣到復(fù)雜邊界區(qū)域;排斥力法一般僅需要在邊界處設(shè)置邊界粒子,避免了在計(jì)算區(qū)域外設(shè)置虛粒子的問題,因此具有很強(qiáng)的適應(yīng)性,但需要構(gòu)造合理的排斥力公式來保證邊界力的計(jì)算精度及穩(wěn)定性;耦合邊界法是排斥力法與虛粒子法的耦合,因此,合理的排斥力模型是耦合邊界法的基礎(chǔ)之一.本文著重研究以下5種典型的排斥力法,流體粒子與邊界粒子的相互作用如圖1所示.

本文采用leap-frog格式的時(shí)間積分方法,即

?表示密度ρ或速度v,xi是粒子i的位置坐標(biāo).為了使計(jì)算過程穩(wěn)定,本文采用Monaghan[7]給出的分別考慮具有黏性耗散和外力作用的時(shí)間步長(zhǎng)表達(dá)式:

圖1 流體粒子與邊界粒子相互作用

方法1

Monaghan[7]最先提出用一組固定在邊界上的粒子對(duì)鄰近的流體粒子施加排斥力,從而防止流體粒子非物理穿透邊界,該排斥力類似于計(jì)算分子力時(shí)的Lennard-Jones方程,表達(dá)式如下:

其中f是作用在單位質(zhì)量上的力,ci代表當(dāng)?shù)芈曀?

2.4 傳統(tǒng)排斥力模型

自由滑移邊界條件的定義如下:

為了提高計(jì)算的精度及穩(wěn)定性,國(guó)內(nèi)外學(xué)者提出了多種方法來施加壁面邊界條件.這些方法大致可分為三類:1)排斥力法;2)虛粒子法;3)耦合邊界法.其中虛粒子法需要將流體粒子映射到邊界外部

其中參數(shù)n1和n2一般取值為12和4,D一般取為最大速度的平方.r0為截止半徑,若其取值過大,會(huì)使鄰近邊界的流體粒子受到較大的邊界力,導(dǎo)致計(jì)算失敗;若取值較小,就很難阻止流體粒子穿透邊界.

方法2

Rogers等[13]提出的排斥力表達(dá)式為

其中nj為邊界粒子 j處的邊界法向,Ψ為流體粒子i到邊界的距離,ξ為xij在粒子 j處切線方向的投影,u⊥為粒子i的速度在nj上的投影.排斥力函數(shù)R(Ψ)的表達(dá)式為

其中?b為兩鄰近邊界粒子的距離,P(ξ)可保證粒子平行于邊界運(yùn)動(dòng)時(shí),排斥力保持不變.ε(z,u⊥)是修正函數(shù),可根據(jù)流體的深度及速度調(diào)整邊界力的大小:

其中

式中h0為初始流體深度,z為流體局部深度,c0為初始聲速.

方法3

為克服Lennard-Jones方程計(jì)算邊界力時(shí)的不足,Monaghan等[15]又提出了一種邊界力公式:

式中K取值為gh0,β一般取4,mi為粒子i的質(zhì)量.

方法4

強(qiáng)洪夫等[14]基于Galerkin形式的SPH方程推導(dǎo)出了新的排斥力公式:

其中Aj是與邊界粒子 j有關(guān)的邊界權(quán)函數(shù),ε為罰參數(shù),其選取是個(gè)較困難的問題,一般需要若干次試算.由于該方法施加的邊界力大小與流體相對(duì)于邊界法向的速度大小有關(guān),當(dāng)流體粒子靠近邊界且相對(duì)速度很小時(shí),會(huì)導(dǎo)致邊界力振蕩,從而導(dǎo)致靠近邊界的流體粒子速度、壓強(qiáng)和密度的不穩(wěn)定.

方法5

與Monaghan提出的邊界力公式相似,Liu等[16]也提出了一種排斥力方法:

上述5種方法各有優(yōu)缺點(diǎn),下文將通過數(shù)值算例進(jìn)行說明.

3 新型排斥力模型

如引言所述,國(guó)內(nèi)外學(xué)者已提出多種排斥力方法,但都存在各自的不足之處.為了克服上述問題,本文基于Galerkin加權(quán)余量法并結(jié)合上述排斥力方法,提出了一種新型排斥力公式(方法6)來施加壁面邊界條件.

Lagrangian描述的動(dòng)量方程的張量形式如下:

vα是速度分量,σαβ是總應(yīng)力張量分量,gα為單位質(zhì)量體積力分量,以上公式中的重復(fù)上標(biāo)符合Einstein求和約定.應(yīng)用Galerkin加權(quán)余量法,動(dòng)量方程(27)式可寫為如下的弱形式:

式中右側(cè)第一項(xiàng)為內(nèi)力項(xiàng),第二項(xiàng)為表面力項(xiàng)即邊界力項(xiàng),第三項(xiàng)為體力項(xiàng),其中?代表整個(gè)計(jì)算域空間,Γ為計(jì)算域邊界區(qū)域,nα為邊界的外法線方向分量,δα代表權(quán)函數(shù)(test function)分量,代表試函數(shù)(trial function)分量.

Galerkin方法的權(quán)函數(shù)及其梯度的表達(dá)式如下:

其中δ(x?xj)為Diracδ函數(shù).

方程(28)右側(cè)的第二項(xiàng)表示邊界力對(duì)動(dòng)量方程的貢獻(xiàn),應(yīng)用罰函數(shù)法對(duì)該項(xiàng)進(jìn)行變換.首先根據(jù)滑移接觸條件(14)來定義邊界的約束罰勢(shì)ΠB,其張量形式的表達(dá)式為

式中ε為罰參量,罰勢(shì)的一階變分形式為

將該罰勢(shì)的變分形式引入到動(dòng)量方程的弱形式(28)式中,并代替邊界作用項(xiàng),得到如下的形式:

將權(quán)函數(shù)和試函數(shù)的近似表達(dá)式(29),(30)和(31)代入(34)式,得到如下形式:

由于δvα的任意性,因此對(duì)于每個(gè)粒子必須滿足下式:

應(yīng)用核函數(shù)的性質(zhì)并對(duì)上式進(jìn)行點(diǎn)積分變化得到如下表達(dá)式:

通過將(39)式與(24)式對(duì)比可確定罰參數(shù)ε及Hij的取值.為使(39)式計(jì)算所得的邊界力與(24)式具有相同的數(shù)量級(jí),本文取

當(dāng)流體靠近壁面時(shí),(39)式所得的邊界力與(vi?j)·nj成正比,為了避免方法4產(chǎn)生的壓力振蕩問題,本文對(duì)該項(xiàng)也進(jìn)行了修正,最終得到的邊界排斥力的定義如下:

4 數(shù)值算例

本節(jié)通過4個(gè)數(shù)值算例來驗(yàn)證本文方法的有效性,并與其他方法進(jìn)行了對(duì)比.為了保證計(jì)算結(jié)果的可對(duì)比性,每種方法都采用相同的粒子配置,核函數(shù)及其影響域也保持一致.

分別通過靜止液柱算例及液柱坍塌算例,對(duì)比了傳統(tǒng)排斥力公式及本文新型排斥力公式的優(yōu)缺點(diǎn),證明了本文方法的有效性.

對(duì)容器中液體靜止算例及潰壩算例,當(dāng)固壁粒子僅采用邊界上的單層粒子分布時(shí),會(huì)造成邊界附近粒子缺失,從而導(dǎo)致邊界附近的流體粒子計(jì)算精度降低,本文采用Liu等[16]提出的虛粒子施加方法來解決該問題.即在邊界粒子外側(cè)設(shè)置兩排虛粒子,邊界粒子的信息通過其鄰近的流體粒子插值得到,虛粒子的信息通過邊界粒子和流體粒子的信息插值獲得,如(41)—(43)式:

通過液體靜止算例,驗(yàn)證了本文方法所得的邊界力與真實(shí)值的一致性;通過潰壩算例,驗(yàn)證了本文方法處理復(fù)雜邊界問題的有效性.

4.1 靜止液柱算例

為了驗(yàn)證邊界力對(duì)流體粒子初始分布的擾動(dòng)特性,本節(jié)選取靜止液柱算例進(jìn)行數(shù)值分析.

該算例中,液柱不受重力作用.其尺寸為0.4 m×0.6 m,邊界長(zhǎng)度為1 m.核函數(shù)統(tǒng)一采用三次樣條函數(shù),粒子i的影響域hi=1.0?di,?di為粒子i的初始尺寸.初始密度ρ0=1000.0 kg/m3,邊界粒子與流體粒子尺寸大小相同,取為0.02 m,時(shí)間步長(zhǎng)為1×10?5s,初始粒子分布如圖2,最底部的一層粒子為邊界粒子,其余粒子為流體粒子.由于流體靜止且不受外力作用,要保證系統(tǒng)動(dòng)量守恒,在計(jì)算過程中流體粒子應(yīng)始終保持靜止.

t=0.3 s時(shí),采用不同方法模擬的結(jié)果如圖3,圖中箭頭方向?yàn)榱Ｗ拥倪\(yùn)動(dòng)方向.從圖3可知,方法4及本文方法能保證流體粒子始終靜止.

圖2 靜止液柱算例初始SPH粒子分布(單位為m)

圖3 t=0.3 s時(shí)用不同方法模擬靜止液柱獲得的SPH粒子分布及速度分布(單位為m) (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

方法1的模擬結(jié)果和截止半徑r0及影響域h的取值有關(guān).若r0＜h/2,則流體粒子保持靜止,否則流體粒子將受到邊界力的擾動(dòng).本文取r0=1.1?d＞h/2,鄰近邊界的流體粒子由于受到邊界粒子排斥力的作用而產(chǎn)生運(yùn)動(dòng);若增大r0的取值,流體粒子所受的排斥力增加,流體粒子的運(yùn)動(dòng)速度也會(huì)增大.選取四次樣條核函數(shù)及高斯核函數(shù)也存在同樣的問題.因此,截止半徑r0的取值問題不利于方法1在實(shí)際問題中的應(yīng)用.

方法2、方法3及方法5中,無論選取何種核函數(shù),與邊界粒子鄰近的流體粒子受到的邊界排斥力均不為0,邊界力的作用使流體粒子產(chǎn)生運(yùn)動(dòng),因此不能保證動(dòng)量守恒;且方法2中邊界排斥力最大,流體粒子運(yùn)動(dòng)速度明顯大于其他方法所得的流體粒子運(yùn)動(dòng)速度.

方法4的罰參數(shù)取值為1.0.在方法4及方法6中,排斥力大小與流體和壁面之間的相對(duì)速度有關(guān),僅當(dāng)流體靠近壁面時(shí),才受到邊界力的作用.這就使初始靜止流體能始終保持靜止,保證了鄰近邊界的流體粒子分布與實(shí)際相符.

若核函數(shù)統(tǒng)一采用高斯核函數(shù)或四次樣條函數(shù),影響域取初始粒子間距的1.0—1.5倍,計(jì)算結(jié)果與上述結(jié)果相似,所得的結(jié)論保持不變.

4.2 液柱坍塌算例

液柱在重力作用下會(huì)向下運(yùn)動(dòng),排斥力的存在將阻止其穿透壁面,本算例對(duì)各種排斥力阻止流體粒子非物理穿透壁面的特性進(jìn)行研究.

核函數(shù)及其影響域、液柱尺寸、初始設(shè)置與上節(jié)相同,但受重力作用,重力加速度g=9.81 m/s2,沿y軸負(fù)方向,

t=0.3 s時(shí)的計(jì)算結(jié)果如圖4所示.從圖4中可知:當(dāng)邊界粒子和流體粒子尺寸大小相同時(shí),方法1、方法3、方法4及方法5并不能阻止流體粒子穿透邊界,究其原因,是由于這4種方法的排斥力大小與邊界粒子和流體粒子之間的距離成反比.要防止流體粒子非物理穿透邊界,有兩種方法:一種是將邊界粒子尺寸設(shè)為流體粒子的1/2或更小,如文獻(xiàn)[15],但是這樣會(huì)增加邊界粒子數(shù)目,降低計(jì)算效率,尤其是三維問題;另一種是使排斥力大小與流體粒子到邊界粒子的法向距離成反比,如方法2和方法6,這種處理方法可減小邊界力精度對(duì)邊界粒子尺寸的依賴,提高計(jì)算效率.

方法4中,罰參數(shù)ε選取需要經(jīng)過若干次試算才能獲得.若罰參數(shù)過大,流體粒子在向下運(yùn)動(dòng)的過程中會(huì)受到過大的邊界力,從而導(dǎo)致該粒子加速遠(yuǎn)離邊界運(yùn)動(dòng),擾亂整個(gè)流場(chǎng)的粒子秩序,致使模擬失效;若罰參數(shù)過小,邊界排斥力很難阻止流體粒子穿透邊界.本算例經(jīng)過3次試算,分別取罰參數(shù)為0.1,1.0,10.0,最終確定罰參數(shù)取1.0.

4.3 容器中液體靜止算例

為了驗(yàn)證本文方法所得的排斥力與邊界力的真實(shí)值是否一致,選取文獻(xiàn)[13]的容器中液體靜止算例進(jìn)行檢驗(yàn).本節(jié)僅選取方法3、方法4與本文方法進(jìn)行對(duì)比,來驗(yàn)證本文方法的有效性.

圖4 t=0.3 s時(shí),不同方法計(jì)算液柱坍塌得到的SPH粒子分布 (a)方法1;(b)方法2;(c)方法3;(d)方法4;(e)方法5;(f)方法6

容器長(zhǎng)2 m,高1 m,水深H=0.9 m.方法3中,計(jì)算參數(shù)與文獻(xiàn)[15]一致,為防止流體粒子穿透邊界,邊界粒子尺寸設(shè)為流體粒子的1/2;方法4中,罰參數(shù)設(shè)為2.0,流體粒子尺寸設(shè)為0.02,邊界粒子尺寸設(shè)為0.015 m;本文方法中所有粒子的尺寸都設(shè)為0.02 m.3種方法的時(shí)間步長(zhǎng)同取2.0×10?5s,所有流體粒子的初始?jí)簭?qiáng)都設(shè)為0,初始密度ρ0=1000.0 kg/m3,重力加速度g=9.81 m/s2,沿y軸負(fù)方向

在重力作用下,流體粒子向下運(yùn)動(dòng),但邊界力的存在阻礙了其進(jìn)一步運(yùn)動(dòng),這樣流體粒子就會(huì)在重力和邊界力的共同作用下振蕩運(yùn)動(dòng),隨著時(shí)間的增加,流體粒子所受的力會(huì)達(dá)到平衡,容器底部的流體粒子的壓強(qiáng)理論值應(yīng)為ρgH,H為液體深度.

t=1 s時(shí),粒子分布狀態(tài)如圖5.方法3的計(jì)算結(jié)果如圖5(a)所示,可知高度越高,靠近邊界的流體粒子離邊界越遠(yuǎn),且該距離大于初始粒子間距,這就導(dǎo)致左上角及右上角處流體粒子的聚集及粒子秩序的紊亂.方法4及本文方法中,邊界力與粒子間的相對(duì)速度有關(guān),從而克服了方法3存在的問題,鄰近邊界的流體粒子具有更好的分布秩序.

圖5 t=1.0 s時(shí)三種不同方法模擬的SPH粒子分布狀態(tài) (a)方法3;(b)方法4;(c)方法6

圖6 容器底部中央的SPH粒子相對(duì)壓強(qiáng)隨時(shí)間變化過程

以容器底部靠近中央的流體粒子為研究對(duì)象,記錄其壓強(qiáng)隨時(shí)間的變化過程如圖6.圖6中壓強(qiáng)為其相對(duì)值,即p/ρgH,液體趨于靜止時(shí),相對(duì)壓強(qiáng)的理論值為1.0.從圖6可知,方法3及本文方法在t=0.6 s時(shí)壓強(qiáng)都達(dá)到了穩(wěn)定狀態(tài),相對(duì)壓強(qiáng)都趨于1.0.本文方法計(jì)算所得的相對(duì)壓強(qiáng)略小于1.0,主要是由人工黏性產(chǎn)生的耗散所導(dǎo)致的.采用方法4時(shí),計(jì)算所得的相對(duì)壓強(qiáng)持續(xù)振蕩,t=0.6 s時(shí)壓強(qiáng)基本趨于穩(wěn)定,但t=1.02 s時(shí),壓強(qiáng)又開始劇烈振蕩.其原因是:當(dāng)鄰近邊界的流體粒子壓強(qiáng)趨于穩(wěn)定時(shí),其運(yùn)動(dòng)速度趨于0,從(16)式可知,邊界力亦趨于0,流體粒子在壓力作用下會(huì)產(chǎn)生較大的加速度及速度向邊界靠近,此時(shí)在(16)式作用下又會(huì)產(chǎn)生較大的邊界力使流體粒子遠(yuǎn)離邊界,這樣鄰近邊界的流體粒子就會(huì)產(chǎn)生振蕩的運(yùn)動(dòng)速度,由(5)式及(4)式可知,粒子運(yùn)動(dòng)速度的劇烈變化會(huì)導(dǎo)致粒子密度及壓強(qiáng)的大幅振蕩.本文方法則克服了方法4存在的缺陷.

圖7 采用本文方法得到的不同時(shí)刻的相對(duì)壓強(qiáng)分布 (a)t=0.8s;(b)t=1.0s;(c)t=1.2s

為了考察壓力場(chǎng)的空間分布特性,圖7給出了0.8,1.0,1.2s三個(gè)不同時(shí)刻的相對(duì)壓強(qiáng)分布狀態(tài).由圖7可知,本文方法計(jì)算所得的壓力場(chǎng)比較光滑;容器底部壓力場(chǎng)穩(wěn)定后的不同時(shí)刻,壓力場(chǎng)的空間分布一致且不隨時(shí)間變化.

4.4 潰壩算例

潰壩是經(jīng)典的自由表面流動(dòng)問題,常用來驗(yàn)證壁面邊界施加方法的有效性.計(jì)算模型如圖8,水柱位于容器左側(cè),右側(cè)有一擋板,初始時(shí)刻水柱靜止,然后將擋板迅速抽出.水柱高度H=2L,寬度為L(zhǎng)=0.146m,容器尺寸為2L×4L.

圖8 潰壩幾何模型

本文對(duì)潰壩過程進(jìn)行了模擬,并與Koshizuka等[18]的實(shí)驗(yàn)以及方法4、方法5所得的結(jié)果進(jìn)行了對(duì)比.計(jì)算參數(shù)為:初始密度ρ0=1000.0 kg/m3,邊界粒子數(shù)目為805,邊界粒子初始尺寸為1.5mm,流體粒子數(shù)目為5000,流體粒子初始間距為 2.92mm,時(shí)間步長(zhǎng) 1×10?5s,Vmax=方法4的罰參數(shù)經(jīng)試算后取0.5.核函數(shù)統(tǒng)一采用三次樣條函數(shù),粒子i的影響域hi=1.0?di.為克服邊界粒子缺失問題,本算例的3種方法都采用了Liu等[16]提出的虛粒子法來提高邊界附近的計(jì)算精度.

3種邊界施加方法所得的模擬結(jié)果如圖9,第1行為0.2,0.6s時(shí)的實(shí)驗(yàn)結(jié)果;第2行為方法4所得的結(jié)果;第3行為方法4所得的結(jié)果;第4行為采用本文方法所得的模擬結(jié)果.圖9中箭頭方向表示粒子的運(yùn)動(dòng)方向,其長(zhǎng)短表示粒子運(yùn)動(dòng)速度的大小.

從圖9中可知,3種方法所得的自由表面運(yùn)動(dòng)過程與實(shí)驗(yàn)結(jié)果都比較符合.但是,采用方法4施加自由滑移邊界條件時(shí),靠近邊界處部分流體粒子的速度明顯大于其周圍區(qū)域流體粒子的速度;隨著時(shí)間的增加,這種非物理的速度振蕩并不會(huì)消失,且邊界流體粒子速度失真會(huì)導(dǎo)致粒子密度和壓強(qiáng)的不穩(wěn)定,最終導(dǎo)致邊界附近流場(chǎng)的數(shù)值模擬失真.采用方法5及本文方法時(shí),靠近邊界的流體粒子始終保持平行壁面邊界運(yùn)動(dòng),粒子速度場(chǎng)分布光滑且秩序良好,很好地克服了方法4所導(dǎo)致的速度振蕩問題.

5 結(jié)論

本文提出一種新型排斥力模型,該模型不需要確定未知參數(shù),可較好地防止流體粒子非物理穿透壁面,能克服流體相對(duì)壁面速度較小時(shí)速度不穩(wěn)定問題.通過液柱靜止算例、液柱坍塌算例及容器中液體靜止算例對(duì)比并驗(yàn)證了本文方法的動(dòng)量守恒性、防止流體粒子非物理穿透邊界的特性及邊界力的精確度;最后,將本文方法應(yīng)用到潰壩算例,進(jìn)一步驗(yàn)證了該方法適應(yīng)復(fù)雜邊界問題的能力.

本文方法可處理復(fù)雜形狀邊界,不需要確定未知參數(shù),可較好地施加固壁邊界條件.本文方法還得出:邊界排斥力的大小要與流體粒子到邊界的法向距離成反比,而不是與流體粒子到邊界粒子的距離成反比,此時(shí)才能保證邊界力不受邊界粒子尺寸的影響,有效防止流體粒子非物理穿透壁面邊界;排斥力模型中需要考慮流體與壁面的相對(duì)速度對(duì)排斥力的影響,否則會(huì)導(dǎo)致邊界附近流體粒子速度失真.

圖9 三種不同邊界施加方法所得的潰壩速度場(chǎng)及與實(shí)驗(yàn)結(jié)果的對(duì)比 (a)實(shí)驗(yàn)結(jié)果,t=0.2 s;(b)實(shí)驗(yàn)結(jié)果,t=0.6 s;(c)方法4,t=0.2 s;(d)方法4,t=0.6 s;(e)方法5,t=0.6 s;(g)本文方法,t=0.2 s;(h)本文方法,t=0.6 s

[1]Lucy L B 1977 Astron.J.82 1013

[2]Gingold R A,Monaghan J J 1977 Mon.Not.R.Astron.Soc.181 375

[3]Zhang A M 2008 Chin.Phys.B 17 927

[4]Sun Z H,Han R J 2008 Chin.Phys.B 17 3185

[5]Zhang A M,Yao X L 2008 Acta Phys.Sin.57 339(in Chinese)[張阿漫,姚熊亮2008物理學(xué)報(bào)57 339]

[6]Monaghan J J 2005 Rep.Prog.Phys.68 1703

[7]Monaghan JJ 1994 J.Comput.Phys.110 399

[8]Libersky L D,Petscheck A G,Carney T C 1993 J.Comput.Phys.109 67

[9]Randles P W,Libersky L D 1996 Comput.Methods Appl.Mech.Eng.138 375

[10]Gotoh H,Sakai T 1999 Coast.Eng.41 303

[11]Liu G R,Gu Y T 2001 Struct.Eng.Mech.246 29

[12]Liu M B,Chang JZ 2010 Acta Phys.Sin.59 3654(in Chinese)[劉謀斌,常建忠2010物理學(xué)報(bào)59 3654]

[13]Rogers B,Dalrymple R 2007 Adv.Num.Model Simul.Tsun.Wave Runup.10 75

[14]Qiang H F,Han Y W,Wang K P,Gao W R 2011 Eng.Mech.28 245(in Chinese)[強(qiáng)洪夫,韓亞偉,王坤鵬,高巍然2011工程力學(xué)28 245]

[15]Monaghan J J,Kajtar J B 2009 Comput.Phys.Commun.180 1811

[16]Liu M B,Shao J R 2011 Sci.China:Technol.Sci.10 1

[17]Chang J Z,Liu M B,Liu H T 2008 Acta Phys.Sin.57 3954(in Chinese)[常建忠,劉謀斌,劉漢濤2008物理學(xué)報(bào)57 3954]

[18]Koshizuka S,Oka Y 1996 Nucl.Sci.Eng.123 421