繆吉倫,陳景秋,張永祥

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

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

1 SPH方法的基本原理

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

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

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

2 流體SPH數(shù)值方法

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

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

質(zhì)量守恒方程

動量守恒方程

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

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

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

2.2 光滑核函數(shù)的選取

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

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

2.3 邊界處理方法

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

2.4 數(shù)值求解

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

3 SPH方法在自由表面流研究中的應(yīng)用

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

3.1 潰壩問題

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

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

3.2 河冰、海冰模擬方面

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

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

3.3 大壩溢流

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

3.4 流固耦合問題

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

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

3.5 水下爆破

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

3.6 多相流的研究

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

3.7 船舶流體力學(xué)

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

3.8 其他領(lǐng)域

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

4 SPH方法的局限及發(fā)展方向

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

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

[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方法的理論及應(yīng)用[J].水利水電技術(shù),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]王剛,李海,季順迎.基于光滑質(zhì)點流體動力學(xué)的海冰熱力-動力數(shù)值模式[J].水科學(xué)進展,2007,18(4):523-530.

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

[12]李大鳴,劉江川,徐亞男.SPH法在大壩表孔泄流數(shù)值模擬中的應(yīng)用[J].水科學(xué)進展,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方法在兩相流動問題中的典型應(yīng)用[J].湖南大學(xué)學(xué)報:自然科學(xué)版,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]鄭興,段文洋.光滑質(zhì)點流體動力學(xué)(SPH)及其算法特性[J].船舶力學(xué),2008(8):550-559.