張力方，張建民

(四川大學(xué) 水力學(xué)與山區(qū)河流開發(fā)保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，成都 610065)

近三十年來，新興的基于無網(wǎng)格法的數(shù)值模擬方法在計(jì)算流體力學(xué)領(lǐng)域發(fā)展十分迅速，較為具有代表性的有光滑粒子流體動力學(xué)和格子玻爾茲曼方法。對于簡單的自由面流動問題，傳統(tǒng)的網(wǎng)格法數(shù)值模擬一般可以得到較為精確的計(jì)算結(jié)果，但是當(dāng)自由液面發(fā)生劇烈變形時，如翻卷、破碎、液滴飛濺等現(xiàn)象，受限于網(wǎng)格尺寸等因素往往無法精確地捕捉。然而無網(wǎng)格粒子法在處理這些問題時能夠體現(xiàn)其極大的優(yōu)勢。潰壩問題作為具有自由面大變形的經(jīng)典問題之一，許多學(xué)者利用這兩種方法對潰壩流進(jìn)行了數(shù)值模擬來驗(yàn)證其可靠性?？娂獋惖萚1]利用SPH方法對下游不同障礙物挑流與消能進(jìn)行了初步分析。Marrone等[2]采用delta-SPH方法研究了潰壩水流沖擊不同形狀障礙物的自由面變形，并將預(yù)測的沖擊壓力與試驗(yàn)值和解析解進(jìn)行了比較，驗(yàn)證了該數(shù)值方法的準(zhǔn)確性和穩(wěn)定性。葉永等[3]利用基于LBM方法研發(fā)的X Flow軟件模擬了有矩形障礙物的三維潰壩問題，結(jié)果表明增強(qiáng)壁面函數(shù)法可以提高計(jì)算的穩(wěn)定性。本文分別將SPH法和LBM法應(yīng)用于模擬二維潰壩、三維潰壩水流碰撞固定障礙物，與試驗(yàn)值對比監(jiān)測點(diǎn)的水深、壓力隨時間變化。而對于潰壩水流沖擊浮體的模擬較少，尤其是LBM方法，因此兩種方法還被應(yīng)用于模擬浮體的位置變化，并分別與試驗(yàn)值對比，分析兩種方法在結(jié)果上的準(zhǔn)確性和異同點(diǎn)，最后對比了計(jì)算效率上的差異，得出兩者的優(yōu)勢和不足，便于在解決實(shí)際問題時選擇合適的方法。

1 研究方法

1.1 光滑粒子流體動力學(xué)

Dualsphysics[4]是基于光滑粒子流體動力學(xué)的一套開源代碼，同時具有CPU并行和GPU并行，可操作性強(qiáng)，能夠與多種算法耦合，目前發(fā)展已經(jīng)較為成熟。該代碼基于WCSPH模型，將流體視為弱可壓縮的，控制方程為拉格朗日形式的連續(xù)性方程和動量方程：

(1)

(2)

式中：t表示時間；u表示速度矢量；ρ表示流體密度；P表示壓力；g為重力加速度；υ為運(yùn)動黏滯系數(shù)。

經(jīng)過核函數(shù)近似和粒子近似將上式表示成SPH形式，應(yīng)用delta-SPH公式[5]，在連續(xù)性方程中引入擴(kuò)散項(xiàng)，并且在動量方程中加入人工黏性項(xiàng)可以消除模擬過程中產(chǎn)生的數(shù)值振蕩，因此控制方程如下所示：

(3)

(4)

式中：ρa(bǔ)、Pa分別為粒子a處的密度、壓強(qiáng)；b表示粒子a的支持域內(nèi)其他粒子；rab表示粒子a和粒子b之間的距離；Wab表示光滑核函數(shù)；δΦ為可調(diào)節(jié)的無量綱參數(shù)，對于大多數(shù)模擬問題推薦取值為0.1。

光滑核函數(shù)和人工黏性項(xiàng)∏ab如下所示：

(5)

光滑核函數(shù)的選擇對于SPH方法的模擬效果至關(guān)重要。核函數(shù)的自變量為相互作用粒子之間的無量綱距離(q)，定義為q=r/h。本文采用的三次樣條函數(shù)，如下所示：

(6)

式中：αD在二維問題中取值為10/7 πh2，在三維問題中取值為1/πh3。

基于弱可壓縮假定，通常使用Tait狀態(tài)方程來確定壓力和密度之間的關(guān)系：

(7)

邊界處理方法采用動力邊界條件，假定邊界粒子滿足與流體粒子相同的控制方程，但是不會根據(jù)作用力進(jìn)行移動，而是保持固定的位置或根據(jù)自定義運(yùn)動函數(shù)進(jìn)行運(yùn)動。當(dāng)流體粒子與邊界粒子之間的距離小于2倍的光滑長度，邊界粒子的密度將增加，從而導(dǎo)致壓力的增加。同時，邊界粒子將對流體粒子產(chǎn)生排斥力，以防出現(xiàn)粒子穿透現(xiàn)象。

1.2 格子玻爾茲曼方法

X Flow軟件基于格子玻爾茲曼方法(LBM)，用介觀模型來模擬流體宏觀行為，將水體離散模擬為微粒，研究其運(yùn)動參數(shù)，具有無需網(wǎng)格、高效并行、邊界條件處理簡單等特點(diǎn)。LBM方法最早是由普通Boltzmann方程發(fā)展而來的。1993年，偽勢模型由Shan和Chen提出[6,7]，對于采用標(biāo)準(zhǔn)的Bhatnagar-Gross-Krook (BGK)碰撞算子的偽勢模型，其相應(yīng)的粒子演化方程為：

fi(x+eiΔt,t+Δt)-fi(x,t)=

(8)

(9)

式中：ρ為宏觀密度；u為宏觀速度；αi為權(quán)重系數(shù)。

2 案例分析

2.1 二維潰壩

如圖1所示，上游水體長L=1 m，高H=2 m，總計(jì)算域大小為長4 m×高3 m，求解尺度為0.01 m，計(jì)算時間為2 s。所有壁面均采用無滑移邊界條件。在SPH中，水體粒子數(shù)為20 000。初始狀態(tài)水體由擋板阻隔，處于靜止?fàn)顟B(tài)，開始計(jì)算時，擋板被瞬間提起，水體在自身重力的作用下坍塌變形。

圖1 二維潰壩模型立視圖(單位：m)Fig.1 Elevation view of the two-dimensional dam-break model

圖2展示了潰壩水流前鋒位置的SPH法和LBM法的計(jì)算結(jié)果，并與Koshizuka[8]試驗(yàn)結(jié)果進(jìn)行對比，圖2中結(jié)果已進(jìn)行無量綱化處理。整體上看數(shù)模結(jié)果與試驗(yàn)相吻合，但隨著水流的行進(jìn)，兩者差距越來越明顯，這是由于在數(shù)值模擬中未考慮底部阻力的影響，計(jì)算出的水流前緣稍快于試驗(yàn)結(jié)果，且SPH方法模擬結(jié)果略優(yōu)于LBM方法。

圖3展示了各時刻潰壩水流的自由面形狀數(shù)值模擬結(jié)果，即使經(jīng)歷了坍塌、碰撞、翻卷、入水等過程，兩種方法均能捕捉到自由面的劇烈變形情況。在t=0.5 s時刻，對于潰壩水流前緣的模擬LBM結(jié)果顯得比較凌亂，而SPH結(jié)果的自由表面相對光滑，但由于SPH的固壁邊界處理不夠完善，在邊界處產(chǎn)生極大的數(shù)值黏性使得粒子黏附在左側(cè)壁面上，自由面發(fā)生滯留；在t=1 s時刻，水體碰撞墻壁向上爬升，SPH結(jié)果的水面線呈一條斜線，水體厚度隨著爬升高度的增加而變薄，LBM結(jié)果的前緣粒子在爬升時逐漸破碎、分離，自由面變得不連續(xù)；在t=1.5 s時刻，中部水體由于重力原因也開始翻卷，SPH結(jié)果中前緣水體由于受到上邊界和下方水體的擠壓，飛濺至較遠(yuǎn)處，而LBM的前緣水體一部分沿著頂部繼續(xù)滑移；在t=2 s時刻，翻卷水體與上游水面融合，兩者結(jié)果均在入水時形成空腔，在劇烈碰撞后，LBM結(jié)果中部分粒子仍散落在空中。綜上所述，兩種方法的模擬結(jié)果與試驗(yàn)現(xiàn)象較為符合，且SPH法非常適合于處理劇烈變形的自由面問題，LBM法在自由面捕捉上仍需稍加改進(jìn)。

圖2 潰壩水流波前位置對比Fig.2 Comparison of wave-front position of dam-break flow

圖3 各時刻下自由液面形狀Fig.3 Time sequence of the free surface shape

2.2 潰壩水流碰撞固定障礙物

在實(shí)際工程中，往往由于復(fù)雜的地形特征需要考慮三維特性，在此開展了三維潰壩水體與下游方形結(jié)構(gòu)物碰撞模擬，如圖4所示?？傆?jì)算域?yàn)?.22 m×1.0 m×1.0 m，上游水體為1.228 m×1.0 m×0.55 m。求解尺度為0.01 m，計(jì)算時間為6 s，其余設(shè)置與上述案例相同。在SPH中，水體粒子數(shù)為669 735。根據(jù)Kleefsman等人[10]在荷蘭海事研究所(MARIN)進(jìn)行的試驗(yàn)，有4個水深測點(diǎn)分布在水池中軸線，一個設(shè)在上游水庫，另外3個設(shè)置在下游空地。在障礙物的前部和頂部上各布置4個壓力傳感器，如圖5所示。

圖4 潰壩水流碰撞固定障礙物模型布置圖(單位：m)Fig.4 Sketch of the dam-break flow against fixed obstacle

圖5 障礙物尺寸及壓力測點(diǎn)布置(單位：m)Fig.5 Sketch of the obstacle and pressure measuring point

通過觀察圖6可知，在SPH 模擬結(jié)果中，t=1 s時刻水流與障礙物碰撞后，自由面產(chǎn)生劇烈變形，兩側(cè)水流繞開障礙物并向中心線運(yùn)動，直至碰撞右側(cè)墻壁，合并形成一個反向回流，與越過障礙物的水流相遇。水流沿X軸正向流速逐漸增大，在右側(cè)墻壁中心線處流速最大達(dá)到3.7 m/s，因?yàn)檫@部分水流未受到阻礙，重力勢能幾乎全部轉(zhuǎn)化為動能。

圖6 在t=1 s時刻SPH方法模擬下自由面變形及速度云圖Fig.6 Free surface deformation and velocity contour simulated by SPH method at t=1 s

圖7展示了H1、H2、H3、H4測點(diǎn)處水深值時歷曲線試驗(yàn)值與模擬值對比，數(shù)值模擬整體趨勢與試驗(yàn)值一致。H1測點(diǎn)離上游最遠(yuǎn)且受障礙物阻擋，水深在約t=1 s后開始迅速增加，水流撞擊右側(cè)墻壁發(fā)生較大變形，試驗(yàn)與數(shù)模結(jié)果均有一定程度的波動，水流重新流向上游，水深減少，當(dāng)上游水流在t=5 s時刻再次流動至H1測點(diǎn)，水深有所回升；H2、H3測點(diǎn)位于障礙物前方，水深值先后開始增加且受障礙物回流影響在短時間內(nèi)陡增，約在t=4.25 s時，來自水庫的回流再次撞擊障礙物，水深值第二次陡增。H4測點(diǎn)位于上游水庫，水深在t=2.5 s降至最低值0.1 m，隨后受回流影響有兩次陡增，一次是來自于撞擊障礙物的回流，另一次是撞擊最右側(cè)墻壁的回流，說明遭遇障礙物后產(chǎn)生的回流會引起水深的不連續(xù)突變，不同于潰壩波的漸變式連續(xù)變化。

圖7 H1、H2、H3、H4測點(diǎn)處水深時歷曲線對比Fig.7 Comparisons of water depth time histories at H1, H2, H3 and H4

圖8展示了P1、P2、P3、P5測點(diǎn)處壓力值時歷曲線試驗(yàn)值與模擬值對比，可以看出數(shù)值模擬較好地捕捉了障礙物受到水流沖擊的壓力變化過程。在t=0.5 s時刻，水流前緣碰撞障礙物前部，壓力急劇上升，隨后緩慢下降，在t=4.5 s時刻壓力曲線出現(xiàn)第二個波峰，是受到反向回流沖擊的影響。SPH由于采用動力邊界條件，邊界粒子與流體粒子之間存在一定間隔，所以測點(diǎn)必須放置在距離邊界1.5倍光滑長度h處，在t=0.5～1.5 s時段出現(xiàn)明顯的壓力曲線波動，并且過高地估計(jì)了壓力峰值，在t=1.5～4.5 s時段壓力曲線比較穩(wěn)定且與試驗(yàn)值一致，對第二個波峰的時間預(yù)測略微延遲，但壓力大小較為吻合。測點(diǎn)P5位于障礙物頂部，約在t=1.5時刻出現(xiàn)壓力峰值，是流體粒子碰撞后一部分翻卷落入上游，另一部分翻越障礙物砸向障礙物頂部，SPH預(yù)測值偏低，在t=4～6 s時刻，出現(xiàn)不合理的負(fù)壓現(xiàn)象，一方面受制于邊界處理方法的缺陷，另一方面在障礙物頂部的水深僅有0.04～0.05 m，而1.5h為0.026 m，邊界處粒子數(shù)偏少，不便于壓力的捕捉。LBM結(jié)果壓力預(yù)測值與試驗(yàn)值吻合良好，無負(fù)壓現(xiàn)象發(fā)生。綜上，對于壓力的處理LBM方法預(yù)測結(jié)果良好。

圖8 P1、P2、P3、P5測點(diǎn)處壓力值時歷曲線對比Fig.8 Comparisons of pressure time histories at P1, P2, P3 and P5

圖9 潰壩水流碰撞可運(yùn)動障礙物模型平面布置圖(單位：m)Fig.9 Plan view of the dam-break flow against transportable obstacle

2.3 潰壩水流碰撞可運(yùn)動障礙物

Amicarelli等人[11]在巴西利卡塔大學(xué)水力學(xué)實(shí)驗(yàn)室進(jìn)行了三維潰壩水流碰撞可運(yùn)動障礙物的試驗(yàn)研究， 模型布置如圖9所示。水槽總長3.5 m，寬0.5 m，上游水庫水深為0.35 m，并在水體右側(cè)設(shè)立閘門，其尺寸為0.037 6 m×0.5 m×0.4 m，閘門在t=2 s時開啟，保持以0.11 m/s的速度向上提起。有兩個固定障礙物布置在下游，分別距左邊界2.5 m和2.95 m，距下邊墻0.06 m和0.29 m，尺寸均為0.15 m×0.15 m×0.75 m。放置在水槽底板上的浮體為一邊長為0.054 m的正方體，質(zhì)量為0.073 kg，密度為464 kg/m3，其質(zhì)心距離左邊界2.532 m，距離上邊墻0.187 m。設(shè)置粒子間距為0.01 m，SPH方法中水粒子數(shù)為86 436，強(qiáng)迫運(yùn)動粒子數(shù)為9 800，浮體粒子數(shù)為252，浮體的動摩擦系數(shù)為0.7，恢復(fù)系數(shù)為0.6，其余邊界粒子動摩擦系數(shù)為0.45，恢復(fù)系數(shù)為0.65。LBM方法中設(shè)置浮體運(yùn)動方式為剛體運(yùn)動，不進(jìn)行任何約束，系數(shù)設(shè)置與SPH方法相同。

試驗(yàn)過程中，潰壩水流挾帶浮體撞向第二個障礙物，由于該浮體密度小于水，被水流提升一定高度，最后移動到右岸，在下游緩慢漂浮，試驗(yàn)監(jiān)測的浮體質(zhì)心坐標(biāo)變化與數(shù)值模擬結(jié)果對比如圖10所示。SPH結(jié)果在t=3～4 s時差異較大，浮體起動時間較試驗(yàn)晚了0.5 s，說明預(yù)測的潰壩水流速度偏小，是由于DEM方法考慮了水流與底板的摩擦，之后與試驗(yàn)結(jié)果吻合較好。LBM方法預(yù)測的浮體起動時間與試驗(yàn)值一致，但是過高估計(jì)了浮體在垂直方向的位移，隨水流被提升至接近0.5 m，在t=4 s時開始下落，受到右側(cè)墻壁的反射波影響，最終漂浮到初始位置的上游，說明LBM方法對于剛體運(yùn)動模擬還不夠完善。

圖10 可運(yùn)動障礙物位置對比Fig.10 Comparison of the position of transportable obstacle

2.4 計(jì)算效率對比

數(shù)值模擬方法的計(jì)算效率也是大多數(shù)學(xué)者們較為關(guān)注的一個重要方面。本小節(jié)將對以上3個案例的計(jì)算時間進(jìn)行對比，如表1所示。程序是在CPU為Intel? CoreTMi7-4790(3.6 GHz)，GPU為NVIDIA GeForce GTX 1080Ti的環(huán)境下運(yùn)行的，但XFlow軟件的LBM方法目前還不支持GPU并行，僅使用CPU并行進(jìn)行計(jì)算。在采用相同CPU并行的條件下，LBM計(jì)算效率略高于SPH方法；采用GPU加速后，SPH 的計(jì)算速度加快了15～25倍，計(jì)算效率顯著提升，說明GPU并行非常適合于粒子法，因此非常有必要實(shí)現(xiàn)LBM方法的GPU并行。

表1 計(jì)算時間對比 s

3 結(jié) 論

SPH和LBM方法在潰壩問題中的模擬結(jié)果與試驗(yàn)值大致吻合，驗(yàn)證了計(jì)算模型的可靠性和穩(wěn)定性。對于自由液面的處理，SPH方法具有較強(qiáng)的自適應(yīng)性使得自由面光滑、清晰，但是會出現(xiàn)粒子附壁現(xiàn)象，而LBM結(jié)果較為凌亂、破碎；對于邊界處壓力值的預(yù)測，LBM方法優(yōu)于SPH方法，SPH法的壓力值預(yù)測出現(xiàn)波動，是由于邊界粒子積分截斷導(dǎo)致的，需要進(jìn)一步對邊界條件處理方法進(jìn)行改善；對于剛體運(yùn)動的模擬，利用耦合DEM的SPH方法具有較好的模擬效果，LBM結(jié)果偏差較大，仍需進(jìn)一步改善；在計(jì)算效率上，使用相同CPU并行的條件下，LBM計(jì)算效率高于SPH，若將GPU并行應(yīng)用于拉格朗日粒子法能夠極大地縮短計(jì)算時間，在計(jì)算工程問題時具有極大的應(yīng)用價值。

□