SPH方法在自由表面流體研究中的應用

繆吉倫,陳景秋,張永祥

(1.重慶交通大學西南水運工程研究所,重慶 400016;2.重慶大學資源環境學院,重慶 400030)

近年來,在流體力學計算領域發展起來一種新型的數值計算方法,即光滑粒子流體動力學法(SPH方法)。SPH方法是一種無網格的純Lagrange方法,1977年由Gingold等[1]提出,最初主要用于解決天體物理學中流體質團無邊界情況下的三維空間任意流動的計算問題。SPH方法無需網格,避免了網格生成的麻煩,而且解決了通常拉氏方法中的網格纏結和扭曲以及網格重劃的問題。近年來,國內外學者采用SPH方法在固壁邊界和自由表面流動等方面進行數值模擬計算,取得了一定的研究成果。

1 SPH方法的基本原理

與傳統的基于網格的方法(如FDM法和FEM法)相比,SPH方法的主要優點在于不需要使用任何提前定義的提供結點連接信息的網格,它用一系列任意分布的粒子質點來代表整個連續介質流體并估計相應的偏微分方程。

SPH方法基于以密度、速度、能量等為變量的偏微分方程組,將描述場的函數用核函數逼近近似表達為任意函數和核函數的乘積的積分:

式中:x與x′為計算域內任意2點的坐標;＜f(x)＞表示坐標 x處的核估計值;D為整個求解區域;f(x′)為坐標 x′處的場量值;h為光滑長度;W(xx′,h)為核函數,它有 2個自變量:粒子間距離和光滑長度h。

2 流體SPH數值方法

2.1 Navier-Stokes方程的SPH粒子近似式

參照SPH方法的粒子近似過程,可以將Navier-Stokes方程寫成以下粒子近似方程[2](不考慮黏度):

質量守恒方程

動量守恒方程

式中:ρi,ρj分別為粒子i和粒子j的密度;t為時間;N為粒子i的支持域中粒子的總數;v為粒子速度;mj為粒子j的質量;Wij為粒子j對粒子i產生影響的光滑函數;α,β表示坐標方向;p為壓力;μ為動力黏性系數;ε為剪切應變率;e為單位質量物質的內能。

根據人工壓縮性,可把一般的不可壓縮流體看作可壓縮流體。通過在水的狀態方程引入壓縮率可模擬具有自由表面的流動:

式中:ρ0是粒子的初始密度;γ為系數,對于水取γ=7;B為用于限制密度的最大改變量,一般為初始壓力。

2.2 光滑核函數的選取

核函數除滿足歸一性和δ函數條件外,還必須具有良好的緊支性,目前常用的核函數有高斯核函數、三次樣條函數、五次樣條函數等。應用較廣的B樣條光滑函數如下[2]:

式中:αd為歸一化常數,在一維、二維、三維情況下分別是1/h,15/(7πh2),3/(2πh3)。

2.3 邊界處理方法

將入口邊界、出口邊界看成特殊自由邊界,界面處的粒子壓力賦值為零或一定的外壓力值。將固壁邊界離散成為邊界粒子,并假定邊界粒子對靠近它的流體粒子施加一個大小適當的中心排斥力,以阻止流體粒子穿越固壁邊界,且排斥力只在近距離上起作用。文獻[3]將邊界力描述為式中:A,n,s為可調參數;r0為粒子初始間距;r為當前時刻邊界粒子和流體粒子的間距;r為長度矢量。當 r＞r0時置 f(r)=0。

2.4 數值求解

對轉換后的常微分方程組進行時間積分,一般采用顯式格式,如預測校正方法、四階 Runge-Kutta方法、標準蛙跳法等,還可用中心差分、Lax-Wendroff差分格式等。對于不可壓縮流體自由表面流動情況,為避免顯式計算產生振蕩,有學者提出用分數步長法分2步求解。為保證解的收斂性,時間步長Δt應滿足CFL條件。

3 SPH方法在自由表面流研究中的應用

SPH方法作為一種無網格方法,在處理移動邊界問題尤其是大規模扭曲和自由表面問題時具有明顯優勢。從目前國內外研究成果來看,SPH方法在流體領域主要應用于以下幾個方面。

3.1 潰壩問題

潰壩問題是國內外研究者比較關注的一個課題,目前的研究多采用FVM,FEM,VOF等方法,而SPH方法由于能逼真復演流體飛濺、融合等復雜自由表面現象而被日益關注。

解決潰壩模擬問題時主要用到Cross液流模型,它能使有效黏性力成為連續變量且避免了數值不穩定性問題,可廣泛應用于牛頓流體及非牛頓流體中。Shao等[4]運用預測校正法(分數步長法)求解了不可壓縮流體的SPH方程組,按時間步長不斷迭代逐步得到不同時刻的粒子位置及流速、壓力等物理量,并模擬出流體碰壁以后的飛濺和融合等現象。經比較,數值模擬結果與解析解非常相近,且隨著粒子數的增加收斂速度加快。Colagrossi等[5]也用SPH方法模擬了潰壩過程,著重分析了波浪前端撞擊壁面后回落到表面的過程。當浪端回落至與上層表面構成封閉環狀空間時,內部氣體被壓縮,使得壓力已不同于大氣壓力,兩相密度的突變是算法不穩定的主要來源,他們將SPH方程中的mj/ρi代之以mj/ρj,避開了這一問題。用修正后的公式成功模擬出流體碰壁回落的過程,與固定網格上解得的Navier-Stokes方程結果一致。

3.2 河冰、海冰模擬方面

目前,大多數海冰與湖冰的動力學模擬采用有限差分法,存在嚴重的數值擴散問題且冰體內部溫度場準確性不理想。尤其當大型冰體邊緣移動時,網格單元會迅速扭曲[6]。隨著純Lagrange粒子方法的產生,這些問題可以迎刃而解。Shen[7]首次將SPH方法應用于河冰模擬,采用SPH耦合FEM方法建立二維數值模型,模擬出復雜的幾何和流動條件下河冰的表面阻塞及動態運輸過程,通過兩相界面處的相互作用將冰動力學及水力學分量結合起來,在水力學方程中考慮了冰層中水的排出量,用非線性滲流理論模擬冰層中的出流,使水的質量守恒,對冰的模擬更加準確。用此模型研究了在密蘇里河下游使用攔冰柵以減少密西西比河河冰輸出的可行性。Oger[8]采用SPH方法模擬了風力作用下浮冰在水面上的運動狀態,并通過引入黏塑性模型考慮了浮冰的流變,數值分析結果與解析解相符。

Ricardo[9]應用SPH方法對海冰進行數值模擬,并利用黏塑性海冰液流模型求解冰場運動方程,從而模擬出斷裂海冰流在楔形渠道內的運動。王剛等[10]引入SPH方法對遼東灣的區域性漂移海冰進行了數值模擬,在采用黏塑性海冰本構方程的基礎上考慮了熱力學因素對其厚度、密集度的影響,得到熱力動力模式的控制方程,研究了海冰厚度和速度矢量的分布規律,模擬出符合實際的海冰演化過程。

3.3 大壩溢流

孫凱[11]對Omkareshwar溢流壩的壩面水流特性進行了研究。李大鳴等[12]將SPH數學模型引入拉西瓦水電站表孔泄流中,提出采用補水邊界的方法來滿足庫區恒定水位條件,模擬了表孔溢流的流場變化及粒子運動過程。計算結果表明光滑粒子可用于模擬高速水流中非線性大變形流動問題。

3.4 流固耦合問題

SPH方法在求解流固耦合問題方面得到了初步應用。Edmond等[13]采用SPH和LES耦合方法對近海岸單波運動機理進行了初步研究,結果顯示,數值模擬的單波剖面與試驗數據以及解析解結果幾乎完全一致;他們還采用SPH結合大渦模型模擬近岸孤立波的爬升,分析了規則及不規則波浪沖越近海構筑物的運動過程,研究成果與相應試驗中的觀測結果非常相近。Mutsuda等[14]將SPH方法應用于沿海岸波浪與結構物相互作用的研究中,數值模擬了波浪與剛性、彈性建筑物相互作用中建筑物產生的運動、變形情況,試驗結果證實了SPH方法的有效性和精確性。Gómez-Gesteira等[15]對波浪沖越壩體時的形狀和幅度進行SPH數值模擬,并對沖擊的位置、速度、壓力進行了比較分析,結果與試驗數據吻合很好。Gotoh等[16]研究了波浪與半淹沒的防波堤的相互作用。

塊體滑坡往往引起自由水面的劇烈變化。Monaghan等[17]研究了矩形塊體沿弧形坡面下滑撞擊入水的過程,介紹了剛體與水相互沖擊的室內試驗和數值模擬結果。2006年,Oger[8]在SPH框架中引入了粒子搜索方法及一種可估計流體對固體邊界作用力大小的方法,用于獲得流固耦合問題中自由運動物體的動力響應,并以不同楔形體滑入水體的試驗結果證實了這種方法的適用性,并指出需要更多的試驗來證實這種方法的可行性,此外還應在求解過程中考慮空氣的運動狀態。

3.5 水下爆破

Liu等[2]應用SPH方法模擬了水下爆破現象。通過對板條和球形TNT炸藥的數值模擬算例,應用SPH方法研究了高能炸藥的真實爆炸模型和人工爆炸模型,同時也對沖擊效應的水介質緩沖進行研究。為了適應水下爆炸沖擊模擬的需要,提出了處理物質交界面的算法,即在那些彼此之間將要發生穿透的粒子上施加懲罰力,并且允許在不同材料和介質的粒子間進行核近似和粒子近似。因此,這種交界面處理方法不僅能通過核近似和粒子近似施加交界面邊界條件,而且能夠有效地防止不同介質粒子之間的非物理穿透。

3.6 多相流的研究

在多相流的研究方面,Monaghan等[18-19]研究了火山噴發中氣塵兩相運動;Cleary[20]研究了多物質熱傳導問題;Richie等[21]研究了多物質星云的形成過程;Colagrossi等[22]、韓旭等[23]在理想氣體狀態方程中引入van der Waals修正項并對粒子運動速度進行修正,研究了密度差異較大時兩相流的流動問題等。近年來,國外有學者開始應用SPH方法研究泥石流、火山流動等問題。

3.7 船舶流體力學

SPH方法用于船舶流體力學僅見于近期的研究。Gómez-Gesteira[24]利用SPH方法對甲板上浪問題進行研究,得到了甲板上浪的水面形態、甲板上流體的流速變化以及甲板的受力狀況;Landrini等[25]用BEM方法結合SPH方法計算船體阻力,船體表面的速度分布采用BEM方法計算,船艏與船艉波的破碎演化過程則采用SPH方法模擬。Benedict等[26]用SPH方法模擬求解有限水深的水翼繞流問題。Souto等[27]模擬了減搖水艙運動與液體所產生的力矩之間的相位滯后問題,對液艙晃蕩所產生的力矩幅值進行了計算,結果表明該方法能夠達到較高的精度。通過這些研究可以看出SPH方法對自由表面波的翻卷、破碎、重入等細節的模擬可以達到較高的精度。在模擬液艙晃蕩問題時,SPH方法計算局部沖擊壓力具有較高的精度,但整體力矩的計算精度不如傳統網格方法的計算精度高。國內還沒有SPH方法用于船舶水動力學的論文報道。

3.8 其他領域

在其他領域方面,Mǘller等[28]通過在 Navier-Stokes方程中引入表面張力項,應用SPH方法模擬了人體組織中的血液流動。Valizadeh等[29]在研究兩相流問題時,對界面粒子質量進行修正,成功模擬了氣泡在流體中的上升及潰滅過程,并復演了2種流體密度差異較大時的潛入現象。Roubtsova等[30]利用SPH方法結合某一地區的三維地形處理洪水演進的范圍,以此確定危險區域。此外,SPH方法也應用于電影特效和游戲中的波浪流動、海嘯對海岸線的沖擊等方面。

4 SPH方法的局限及發展方向

由于SPH方法無需生成網格,避免了大量的單元劃分,克服了有限元方法局部近似所引起的誤差,能夠適應扭曲變形,在大變形問題中自適應性強。SPH方法近年來在流體研究領域發展迅速,但也存在以下局限[31-32]:①光滑函數的選擇尚無明確標準。光滑函數決定了近似模式和粒子的有效支持域,特別是中心峰值對SPH方法的近似精度起著重要的作用。目前已經提出了多個光滑核函數形式,但針對某具體問題應選擇哪種核函數并無明確標準。②邊界粒子的處理方法存在缺陷。邊界上或臨近邊界的粒子在積分時會被邊界截斷,故SPH方法不能完全適用于整個區域。雖然目前采用了虛粒子法改進,邊界粒子結果仍不理想。③計算工作量大,影響運行速度。在計算中可能出現張力不穩定性振蕩,導致計算不收斂。

針對SPH方法的內在缺陷,一些學者提出了相應的修正或改進方法,SPH方法的穩定性、精度、收斂性也提高很多。然而,目前的分析大多基于均勻分布的粒子,對于如粒子高度分布無序、粒子的光滑長度有變化等更普遍的情況尚有待研究。流場邊界附近粒子運動的壓力特性,其計算結果的可靠性尚未驗證。在邊界處理中,如何引入數學表達使邊界的精度得到提高,從而避免由于邊界不穩定導致擴散還需進一步研究,對SPH方法耗時的臨近粒子搜索計算問題需進行進一步的討論,以滿足實時試驗的要求。總之,雖然SPH方法在諸多領域得到了應用,但尚未形成完善的理論體系,目前應鞏固SPH方法的理論基礎,在修正其內在缺陷方面作進一步研究。

[1]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:375-389.

[2]LIU Gui-rong, LIU Mou-bing.Smoothed particle hydrodynamics:a meshfree particle method[M].Singapore:World Scientific Publishing Company,2003:26-32.

[4]SHAO Song-dong,GOTOH H.Pressure analysis of dam-break and wave-breakingby SPH model[J].Annual Journal of Hydraulic Engineering,JSCE,2003,47:403-408.

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

[6]孫曉艷,王軍.SPH方法的理論及應用[J].水利水電技術,2007(3):44-46.

[7]SHEN Hung-tao.SPH simulation of river ice dynamics[J].Journalof Computational Physics,2000,165(2):752-770.

[8]OGER G.Two-dimensional SPH simulation of wedge water entries[J].Journal of Computational Physics,2006,213(2):803-822.

[9]RICARDO G.Flow of fractured ice through wedge shaped channels:smoothed particle hydrodynamics and discrete element simulations[J].Mechanics ofmAterials,1998,29(12):1-17.

[10]王剛,李海,季順迎.基于光滑質點流體動力學的海冰熱力-動力數值模式[J].水科學進展,2007,18(4):523-530.

[11]孫凱.基于SPH方法對自由表面流的數值模擬及可視化研究[D].鄭州:鄭州大學,2007.

[12]李大鳴,劉江川,徐亞男.SPH法在大壩表孔泄流數值模擬中的應用[J].水科學進展,2008,19(6):841-844.

[13]EDMOND Y M L,SHAO Song-dong.Simulation of near-shore solitary wave mechanics by an incompressible SPH method[J].Applied Ocean Reasearch,2002,24:275-286.

[14]MUTSUDA H,SHINKURA Y.An eulerian scheme with lagrangian particlesforsolving impactpressure caused by wave breaking[C]//Proceedings of the Eighteenth International Offshore and Polar Engineering Conference.Vancovver:ISOPE,2008:162-165.

[15]GóMEZ-GESTEIRA M,DALRYMPLE R A.Using a 3D SPH method for wave impact on a tall structure[J].Journal of Waterway Port Coastal and Ocean Engineering,ASCE,2004,130(2):63-79.

[16]GOTOH H,SHAO Song-dong.SPH-LES modelfor numerical investigation of wave interaction with partially immersed breakwater[J].Coast Eng J,2004,46(1):39-63.

[17]MONAGHAN J J,KOS A,ISSA N.Fluid motion generated by impact[J].Journal of Waterway Port Coastal and Ocean Engineering,ASCE,2003,129(6):250-259.

[18]MONAGHAN J J,KOCHARYAN A.SPH simulation ofmultiphase flow[J].Computer Physics Commun,1995,87(1/2):225-235.

[19]MONAGHAN J J.Implicit SPH drag and dusty gas dynamics[J].Journal of Computational Physics,1997,138(2):801-820.

[20]CLEARY P W.Modelling confined multimaterial heat andmAss flows using SPH[J].AppliedmAthematical Modelling,1998,22:981-993.

[21]RICHIE B W,THOMAS P A.Multiphase smoothed particle hydrodynamics[J].Mon Not R Astron Soc,2001,323(3):743-756.

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

[23]韓旭,楊剛,龍述堯.SPH方法在兩相流動問題中的典型應用[J].湖南大學學報:自然科學版,2007,34(1):28-32.

[24]GóMEZ-GESTEIRA M.Green water overtopping analyzed with a SPH model[J].Ocean Engineering,2005,32:223-238.

[25]LANDRINIM,COLAGROSSI A,TULIN M P.Breaking bow and stern waves:numerical simulations[EB/OL].[2001-04-22].http://www.iwwwfb.org/Abstracts/iwwwfb16/iwwwfb16-23.pdf.

[26]BENEDICT D R,DALRYMPLE R A,GESTEIRA M,et al.SPH for naval hydro-dynamics[EB/OL].[2003-04-06].http://www.iwwwfb.org/Abstracts/iwwwfb18/iwwwfb18-43.pdf.

[27]SOUTO A,DELORME L,PéREZ-ROJAS L,et al.Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics[J].Ocean Engineering,2006,33:1462-1484.

[28]MüLLER M,SCHIRM S,TESCHNER M.Interactive blood simulation for virtual surgery based on smoothed particle hydrodynamics[J].Technology and Health Care,2004(12):25-31.

[29]VALIZADEH A,SHAFIEEFAR M,MONAGHAN J.Modeling two-phase flows using SPH method[J].Journal of Applied Sciences,2008,8(21):3817-3826.

[30]ROUBTSOVA V,KAHAVITA R.The SPH technique applied to free surface flows[J].Computers&Fluids,2006,35(10):1359-1371.

[31]SWEDE J W,HICKS D L,ATTAWAY S W.Smoothed particle hydrodynamics stability analysis[J].Journal of Computational Physics,1995,116(1):123-134.

[32]鄭興,段文洋.光滑質點流體動力學(SPH)及其算法特性[J].船舶力學,2008(8):550-559.