弱可壓δ-SPH算法在潰壩流中的應(yīng)用

高仁祖,顧聲龍 ,牛恒崗

(1.青海大學(xué)水利電力學(xué)院，青海 西寧 810016； 2.黃河上游生態(tài)保護(hù)與高質(zhì)量發(fā)展實驗室，青海 西寧 810016)

光滑粒子流體動力學(xué)(Smoothed Particle Hydrodynamics，SPH)是一種基于光滑過程的無網(wǎng)格拉格朗日格式。該方法將計算區(qū)域(材料)離散為一系列的粒子，這些粒子具有質(zhì)量、體積、速度、加速度和溫度等物理量，粒子的集合表達(dá)了材料區(qū)域的物理量，而粒子之間通過核函數(shù)相互作用。SPH方法首次被Gingold等[1]提出，并在流體動力學(xué)中運(yùn)用。近20年來，由于SPH方法在解決復(fù)雜的自由表面流上有獨(dú)特的優(yōu)勢，在潰壩流動[2]、圓柱繞流[3]、近海岸波浪破碎[4]、多相流[5]、流體—結(jié)構(gòu)相互作用[6]、河流冰動力[7]、淺水應(yīng)用方面[8-9]等方面得到廣泛應(yīng)用。

近年來，許多研究人員對SPH方法進(jìn)行完善和改進(jìn)，特別是在計算精度及計算效率方面不斷探索，試圖得到更為準(zhǔn)確的模擬結(jié)果，以便在更多的工程中進(jìn)行應(yīng)用。SPH方法通常分為弱可壓SPH方法和不可壓SPH方法。對于前者，會使用狀態(tài)方程把壓力場和密度場聯(lián)系起來[10]；對于后者，壓力是通過求解Poisson方程來確定的[11]。通常，弱可壓格式更適合自由表面流，因為沿自由表面的邊界條件是隱式滿足的，這意味著在流動變化過程中，弱可壓格式不需要對自由表面進(jìn)行顯式檢測。但是，弱可壓格式的一個主要缺點(diǎn)是在壓力和密度場中會產(chǎn)生虛假的數(shù)值振蕩。為了消除或者減少虛假的數(shù)值振蕩，Monaghan等[12]在一維激波問題中引入了一個基于Neumann-Richtmeyer人工黏性的人工黏性項，為了抑制沖擊波中出現(xiàn)的振蕩，這種人工黏性項有助于穩(wěn)定SPH算法和減少其虛假的振蕩。Andrea等[13]提出了一種人工黏性項修正的形式，用最小二乘法積分插值來濾除密度場，并在程序中給出了結(jié)果。Ferrari等[14]在Rusanov通量的基礎(chǔ)上建立了一個數(shù)值擴(kuò)散項，這一項已被添加到連續(xù)性方程中，以減少密度場內(nèi)部的數(shù)值噪聲。Molteni等[15]也提出了相應(yīng)的數(shù)值擴(kuò)散項，并在連續(xù)性方程中運(yùn)用。這些擴(kuò)散項的加入給出了很好的計算結(jié)果，但是與流體靜力解不符。Molteni等[15]通過引入一個閾值密度來避免這個問題，這樣擴(kuò)散項只在壓力超過流體靜力解時起作用，可是這種策略導(dǎo)致了擴(kuò)散項作用的急劇減少。

基于以上問題，Antuono等[16]對Molteni等[15]提出的擴(kuò)散項進(jìn)行修正，重新定義了一個新的方程組，其擴(kuò)散項都包含在連續(xù)性方程和能量方程中，其黏性項包含在動量方程中，其中擴(kuò)散項系數(shù)δ和黏性項系數(shù)α均是人為可調(diào)的參數(shù)，目的是解決密度不連續(xù)所帶來的數(shù)值計算不穩(wěn)定問題，以及減小虛假的高頻噪聲；其動量方程與具有人工黏性項的標(biāo)準(zhǔn)弱可壓SPH格式一致，被證明與流體靜力解相符，并適當(dāng)?shù)叵藟毫龊兔芏葓龅奶摷贁?shù)值振蕩。

本文基于Antuono等[16]提出的具有擴(kuò)散項的SPH格式，自主開發(fā)δ-SPH算法程序，探討與研究具有該擴(kuò)散項的SPH方法在自由表面流方面的應(yīng)用。

1 模擬方法

1.1 質(zhì)量守恒方程

在拉格朗日型式下的弱可壓牛頓流體的連續(xù)性方程表示為：

(1)

方程(1)中，u為速度矢量，ρ為流體密度，p是關(guān)于密度變化的壓力場。

弱可壓SPH模型中，狀態(tài)方程[17]為：

(2)

方程(2)中，c0和ρ0分別是流體的聲速與初始流體密度。為滿足弱可壓條件，要求密度的波動范圍小于1%ρ0[17]。在流體計算域Ω內(nèi)，需要設(shè)置全局初始聲速c0，并且滿足[18]:

(3)

通過推導(dǎo)將連續(xù)性方程轉(zhuǎn)化為離散化的SPH形式[17]:

(4)

連續(xù)性方程(4)中，擴(kuò)散項Dab可以表示為：

(5)

(6)

(7)

方程(7)中，下標(biāo)a、b代表計算域內(nèi)任意2個粒子，?a表示相對于粒子a坐標(biāo)的梯度算子，Wab表示光滑核函數(shù)(本文使用Wendland C2核函數(shù)[19])，ρa(bǔ)和ua分別是第a個粒子的密度、速度矢量，Va是a粒子影響域的面積(本文計算域Ω是二維的)，δ是人為可調(diào)的擴(kuò)散項系數(shù)。

1.2 動量守恒方程

根據(jù)動量守恒定律，流體的動量方程如下：

(8)

(9)

在方程(8)和方程(9)中，g和r分別代表流體的重力加速度和位置矢量，p是關(guān)于密度變化的壓力場，V是黏性應(yīng)力張量，對于牛頓流體，可以表示為：

V=λtr(D)I+2μD

(10)

其離散化的SPH形式[17]為：

(11)

(12)

(13)

上式中，pa、ga和ra分別是第a個粒子的壓力、重力加速度和位置矢量，α是人為可調(diào)的黏性項系數(shù)。

圖1 程序流程圖Fig.1 Program flow chart

2 模擬工況設(shè)計

2.1 計算流程

開發(fā)的C++δ-SPH模型程序的基本程序流程圖如圖1所示。流程第一步：讀取幾何模型，并對幾何模型(計算域)進(jìn)行粒子化處理；流程第二步：選取或者設(shè)置光滑核函數(shù)(核函數(shù)的選取只是在程序里選擇哪一種核函數(shù)來參與計算，并沒有提前遍歷求解出每個粒子的核函數(shù)及偏導(dǎo)數(shù)，而是每次遍歷粒子需要時，才去求解該粒子所需的核函數(shù)及偏導(dǎo)數(shù))；流程第三步：對每個粒子進(jìn)行粒子配對，其鄰近粒子搜索用到了鏈表搜索法；流程第四步：遍歷每個粒子，通過上述公式(6)進(jìn)行密度梯度修正，以確保在自由表面處核函數(shù)粒子缺失條件下，密度梯度具有二階精度，此步的計算結(jié)果參與第六步的計算；流程第五步：通過虛粒子法進(jìn)行固壁邊界處理，其中遍歷每個虛粒子，通過狀態(tài)方程反推出它們的密度；流程第六步：遍歷每個流體粒子，求解上述方程(4)和方程(11)的左項；流程第七步：時間進(jìn)程采用蛙跳法推進(jìn)每一步求解，并且每500步輸出1次文件；流程最終需要進(jìn)行計算時間的判斷，直到計算達(dá)到要求。

2.2 邊界條件的處理

SPH方法中，靠近邊界的流體粒子支持域內(nèi)只有部分粒子，在固壁邊界外無流體粒子，出現(xiàn)了積分被邊界截斷的問題，造成邊界附近的粒子穿過邊界，因此需對固壁邊界進(jìn)行處理，以獲得正確的計算結(jié)果。本文使用Adami等[20]提出的一種在流體外布置多層虛粒子實施固壁邊界的方法，這種方法直接將內(nèi)部流體粒子的壓力插值到虛粒子上，如方程(14)，再利用狀態(tài)方程反推出虛粒子的密度，虛粒子通過壓力梯度對流體施加邊界力，該邊界處理方法對復(fù)雜幾何邊界問題應(yīng)用很好。

(14)

上式中，ab表示虛粒子的加速度。

2.3 經(jīng)典潰壩案例

本研究采用經(jīng)典潰壩案例進(jìn)行測試，并與Lobovsky等[21]研究得出的實驗值做對比分析，數(shù)值模擬與實驗的比尺Lb=3。Lobovsky等[21]的實驗?zāi)Ｐ鸵运淅锩娴乃w布置設(shè)置了2組工況，分別為水深600 mm和300 mm，其中水箱長度設(shè)置為1 610 mm，高為600 mm，水體長為600 mm。本文選取工況為水深300 mm 進(jìn)行對比分析，數(shù)值模型初始布置如圖2所示。初始的水箱長度設(shè)置為L=537 mm，在水箱里布置長L1為200 mm和d為100 mm的矩形流體區(qū)域。粒子初始間距為1 mm，流體粒子總數(shù)為2萬個，時間步長為10-4s，模擬總時長為10 s。圖2中橫、縱坐標(biāo)以量綱化X/d、Z/d的形式給出。圖2中右下角P1和P2點(diǎn)為壓力監(jiān)測點(diǎn)，其坐標(biāo)分別為5.42,0.05和5.42,0.32；在Lobovsky等[21]實驗中對應(yīng)的監(jiān)測點(diǎn)為Sensor 1和Sensor 4。

圖2 潰壩數(shù)值模擬模型示意圖Fig.2 Schematic diagram of numerical simulation model of dam break

圖3 數(shù)值模擬模型示意圖Fig.3 Schematic diagram of numerical simulation model

2.4 帶有三角形障礙物的經(jīng)典潰壩案例

由于沒有實際潰壩工況，本研究設(shè)計了帶有三角形障礙物的經(jīng)典潰壩案例，在標(biāo)準(zhǔn)SPH中加入擴(kuò)散項之后進(jìn)行測試。數(shù)值模型的初始布置如圖3所示。初始的水箱長度設(shè)置為L=2.05，在水箱里布置長L1為0.5和d為1.0的矩形流體區(qū)域。粒子初始間距為0.01，流體粒子總數(shù)為5 000個，時間步長為10-4s，模擬總時長為10 s。

3 結(jié)果及分析

3.1 經(jīng)典潰壩案例模擬結(jié)果的分析

在此次測試案例中，選取了黏性項系數(shù)α為0.02，并考慮標(biāo)準(zhǔn)SPH格式(即沒有擴(kuò)散項的SPH模型)與具有擴(kuò)散項的SPH格式(擴(kuò)散項系數(shù)δ取0.1)。為了研究壓力場和密度場的虛假數(shù)值振蕩問題，在狀態(tài)方程中，壓強(qiáng)是密度的函數(shù)，因此，本研究數(shù)值模擬的參數(shù)以壓強(qiáng)為主，最終通過壓強(qiáng)的波動變化來分析其虛假數(shù)值振蕩問題。圖4給出了兩種SPH格式下不同時間的壓強(qiáng)變化云圖，時間以無量綱的形式表示為t(g/d)1/2，監(jiān)測點(diǎn)處的壓強(qiáng)以無量綱的形式表示為P/ρgd。從圖4a組與圖4b組的對比發(fā)現(xiàn)，在圖4a組左圖的近壁面流體壓強(qiáng)色值及右圖的壓強(qiáng)變化色值明顯比較紊亂，其主要原因是在標(biāo)準(zhǔn)SPH格式下，數(shù)值在模擬的過程中產(chǎn)生了虛假的數(shù)值振蕩，為此在連續(xù)性方程中考慮Antuono等[16]的擴(kuò)散項，得到圖4b組壓強(qiáng)變化云圖。在兩種不同時間下給出的壓強(qiáng)變化云圖結(jié)果明顯改善，壓強(qiáng)色值變化明顯，比圖4a組呈現(xiàn)規(guī)律分布。

圖4 t(g/d)1/2=2.81與t(g/d)1/2=6.59時的壓強(qiáng)變化云圖Fig.4 Nephogram of pressure variation when t(g/d)1/2=2.81 and t(g/d)1/2=6.59

圖5 監(jiān)測點(diǎn)P2處隨時間變化的壓強(qiáng)曲線分布圖Fig.5 Pressure curve of monitoring point P2 with time

圖5給出了監(jiān)測點(diǎn)P2處隨時間變化的壓強(qiáng)曲線分布圖。紅點(diǎn)曲線表示Lobovsky等[21]的實驗值；黑實線代表原有的標(biāo)準(zhǔn)SPH格式數(shù)值模擬得到的壓強(qiáng)檢測結(jié)果，時間t(g/d)1/2范圍為3.0～6.0時，數(shù)值上下振蕩較為嚴(yán)重，監(jiān)測點(diǎn)處的壓強(qiáng)振蕩范圍為0.3～1.3；綠實線代表經(jīng)過擴(kuò)散項處理后的SPH數(shù)值計算的結(jié)果。在原有弱可壓SPH格式的連續(xù)性方程中加入擴(kuò)散項后，監(jiān)測點(diǎn)處的壓強(qiáng)振蕩范圍為0.4～0.8，明顯削弱了虛假的數(shù)值振蕩，波動范圍的削減意味著降低了弱可壓SPH中的密度和壓力場虛假數(shù)值振蕩的影響。

圖6 監(jiān)測點(diǎn)P1處隨時間變化的壓強(qiáng)曲線分布圖Fig.6 Pressure curve of monitoring point P1 with time

圖6為監(jiān)測點(diǎn)P1處隨時間變化的壓強(qiáng)曲線分布圖。對基于δ-SPH方法得到的數(shù)據(jù)進(jìn)行了曲線擬合處理，從圖6發(fā)現(xiàn)，具有擴(kuò)散項的SPH格式所得到的數(shù)值模擬結(jié)果與Lobovsky等[21]的實驗值吻合較好，在t(g/d)1/2=2.4時，壓強(qiáng)曲線呈突增趨勢，隨后又緩慢下降，此時，流體流動處于潰壩后碰撞另一側(cè)壁面及沿壁面向上流動的狀態(tài)。數(shù)值模擬結(jié)果與實驗結(jié)果的平均誤差為7%。本實驗采用誤差分析標(biāo)準(zhǔn)，其誤差分析是根據(jù)監(jiān)測點(diǎn)處各段時刻取實驗值及模擬值作為誤差，最后對這些誤差進(jìn)行求和、求平均得到最終的平均誤差。

綜上所述，在連續(xù)性方程中考慮Antuono等[16]的擴(kuò)散項之后，有效地改善了數(shù)值模擬結(jié)果。在壓力監(jiān)測點(diǎn)的曲線圖中，實驗值變化趨勢和數(shù)值變化趨勢一致，說明本次開發(fā)的算法程序可以很好地應(yīng)用于自由表面流方面。

3.2 帶有三角形障礙物的經(jīng)典潰壩案例模擬結(jié)果的分析

圖7給出了t=0.26 s與t=0.46 s時的壓強(qiáng)變化云圖。圖7a組圖通過標(biāo)準(zhǔn)的SPH格式數(shù)值模擬得到。從圖7a中發(fā)現(xiàn)，左圖壓強(qiáng)色值變化分層不明顯，局部色值及右圖延伸出來的水舌部分色值明顯紊亂，局部的壓強(qiáng)大于0值，這些現(xiàn)象都是由虛假的數(shù)值振蕩造成的。圖7b組圖是具有擴(kuò)散項的SPH格式數(shù)值模擬得到。從圖7b可以清晰看出，壓強(qiáng)色值分層明顯，且水舌的色值基本趨于0值，相比于圖7a組圖，在SPH中加入擴(kuò)散項后，壓強(qiáng)云圖效果明顯改善。

圖7 t=0.26 s與t=0.46 s時的壓強(qiáng)變化云圖Fig.7 Nephogram of pressure variation when t=0.26 s and t=0.46 s

4 討論與結(jié)論

近些年，相關(guān)研究不斷地改進(jìn)SPH方法，尤其對于弱可壓SPH方法。但在應(yīng)用于自由表面流方面，連續(xù)性方程帶來了數(shù)值振蕩問題。為此，Monaghan等[12]、Andrea等[13]、Ferrari等[14]、Molteni等[15]、Antuono等[16]對此類數(shù)值振蕩問題進(jìn)行了研究，其中Antuono等[16]在弱可壓SPH方法中解決數(shù)值振蕩問題所給出的新格式δ-SPH最為先進(jìn)，在連續(xù)性方程及能量守恒方程中加入擴(kuò)散項，在動量守恒方程中加入了黏性項，其中擴(kuò)散項系數(shù)δ和黏性項系數(shù)α對于不同流體流動情況可以進(jìn)行人為調(diào)節(jié)，能夠促使格式穩(wěn)定，有效地解決了數(shù)值振蕩問題，與流體靜力解相符。

本文通過開發(fā)δ-SPH算法程序，對經(jīng)典潰壩案例與帶有三角形障礙物的潰壩案例進(jìn)行數(shù)值模擬，得出以下結(jié)論：

(1)通過對δ-SPH算法的數(shù)值模擬結(jié)果與文獻(xiàn)[21]實驗結(jié)果的分析比較表明，標(biāo)準(zhǔn)的SPH格式數(shù)值模擬時受到虛假數(shù)值振蕩的影響，加入人為可調(diào)系數(shù)的擴(kuò)散項之后，數(shù)值模擬的結(jié)果明顯改善。尤其在數(shù)值模擬三角形障礙物的潰壩案例中，擴(kuò)散項的加入對模擬效果的改善更為直觀。

(2)從圖5和圖6中發(fā)現(xiàn)，通過監(jiān)測點(diǎn)P1處的壓強(qiáng)檢測數(shù)據(jù)，對數(shù)值模擬結(jié)果與實驗結(jié)果進(jìn)行對比分析，其變化趨勢一致，平均誤差在7%左右，驗證了該算法具有很好的穩(wěn)定性，體現(xiàn)了弱可壓δ-SPH算法能很好地在潰壩流中進(jìn)行應(yīng)用。

(3)圖5的標(biāo)準(zhǔn)SPH格式和δ-SPH格式數(shù)值模擬的壓強(qiáng)隨時間的波動變化說明了本次開發(fā)的δ-SPH算法，在一定程度上降低了虛假的數(shù)值高頻波動，其改善效果明顯。