基于光滑粒子流體動力學方法的潰壩水流模擬

劉慧玲 李海橋

摘要：受全球氣候變暖的影響，由極端天氣引發(fā)的類似潰壩等問題發(fā)生的概率大大增加，深入研究潰壩水流的水動力特性勢在必行。在分析光滑粒子流體動力學基本原理的基礎上，提出了一種改進的邊界處理方法，即將接近壁面的流體視為層流，在耦合動力邊界附近引入層流黏性近似邊界層理論。采用該方法對潰壩水流進行數(shù)值模擬，將SPH數(shù)值模擬得到的外輪廓、自由液面高度以及壓力與實驗結果進行了比較和分析。結果表明：改進的邊界處理方法較完整地得到了水流與壁面相互作用而產生的多種復雜的物理現(xiàn)象，其外部輪廓與實驗非常吻合;自由表面融合過程中液面間沖擊的能量耗散會導致融合后的液面高度存在一些差異;不同監(jiān)測點處壓力隨時間的變化基本落在置信區(qū)間之內。數(shù)值模擬結果與實驗結果吻合度較高，驗證了改進方案的可靠性和計算結果的準確性。

關?鍵?詞：潰壩水流模擬; 自由液面; 耦合動力邊界; 光滑粒子流體動力學（SPH）

潰壩通常指由自然災害或人為因素引起的蓄水庫圍堰突然坍塌，造成人類生命、生態(tài)環(huán)境和生產經濟巨大損失的災害性事件。隨著氣候變暖的影響，極端天氣出現(xiàn)的概率越來越大，由極端天氣導致的類似于潰壩問題出現(xiàn)的概率也大大增加，給人們的生產生活帶來了極大影響。研究潰壩水流的動態(tài)演變過程有利于理解洪峰水動力傳播特性，從而為沿岸防災減災策略的制定提供科學依據(jù)。因此，對潰壩問題的探索具有顯著的科學研究意義和工程應用價值[1]。

目前對潰壩問題的研究主要有3種方法：理論分析、數(shù)值模擬和模型實驗[1-2]。早期對潰壩問題的研究主要以理論分析為主，但理論研究通常局限于比較簡單的情形。近年來，隨著測試技術的不斷更新，潰壩問題的實驗研究也從簡單逐步趨于復雜[2-4]，Lobovsky等對潰壩實驗做了100次的重復實驗[2]，由于該實驗數(shù)據(jù)詳實，受到學術界的高度關注[5-7]。但是由于實驗初始條件存在一定的隨機性，100組實驗數(shù)據(jù)僅展示了液面高度、監(jiān)測點壓力等重要動態(tài)參數(shù)的變化趨勢[8-10]，也未能給出經驗或半經驗公式[11]，這也說明目前的實驗操作水平依然存在很多局限性。近些年來，隨著計算機技術的飛速發(fā)展和偏微分方程求解方法的改進，諸多數(shù)值計算方法廣泛地應用于對潰壩問題的數(shù)值研究中[2-4，12-14]，并逐步將潰壩問題作為一個經典的數(shù)值驗證算例。

由于潰壩問題伴隨著復雜的自由表面演化和流體大變形的強非線性流體運動過程，利用拉格朗日型網格方法，如有限元法（FEM）等處理流體域的大變形時依然會遭受網格的畸變和纏結等因素影響，使計算精度降低和計算成本增加;歐拉型網格方法，如有限差分法（FDM）、有限體積法（FVM）在追蹤復雜自由表面的演化過程時還存在一定的限制[15-16]，雖然VOF[6-7]，Level-set[17]，尤其是耦合VOF和Level-set的方法能較好地模擬該類問題[18]，但對網格精細程度要求依然很高。而無網格法，尤其是拉格朗日型的無網格方法恰好避免了上述問題，在涉及復雜自由表面演化的非線性流動預測中應用越來越廣泛，其中最典型的是光滑粒子流體動力學（Smoothed Particle Hydrodynamics，SPH）方法[15-16，19-20]。SPH方法作為一種純粹的粒子方法，具有完整的拉格朗日粒子性質，在處理伴隨自由表面演化過程的大變形問題中，既能表征自由表面屬性，又能準確描述流體的輸運過程，表現(xiàn)出獨特的優(yōu)勢。

在以往的潰壩問題數(shù)值模擬過程中，固壁邊界的數(shù)值處理方法對計算精度有非常重要的影響[13]。而在真實的流體流動過程中，在固壁邊界附近有明顯的邊界層，邊界層的厚度與主流區(qū)域的紊亂程度相關，但在邊界層區(qū)域均表現(xiàn)出層流流動的特征。目前的潰壩流動實驗研究中，實驗裝置均是規(guī)則的玻璃容器，表面較光滑，且潰壩水流還沒有發(fā)展到整個流體區(qū)域為湍流的程度，根據(jù)邊界層理論，此時的邊界層主要為層流邊界層。本文在計算過程中，將接近壁面的流體視為層流，加入層流黏性來近似概化黏性底層和邊界層，計算潰壩流動，并與實驗研究進行對比。

1?光滑粒子流體動力學

1.1?基本原理

光滑粒子流體動力學（SPH）是一種純拉格朗日型無網格粒子方法，它使用一組粒子離散和代表所模擬的介質（流體或固體），并且基于粒子體系近似和計算介質運動的控制方程。然而不同于一般無網格方法，SPH方法中粒子不僅用于計算場變量與近似控制方程，同時也代表介質系統(tǒng)，具有諸如密度、壓力、速度、內能等宏觀物理量，相當于物質點，并且能夠伴隨介質的運動以當?shù)厮俣纫苿?。根?jù)SPH方法的基本理論，其控制方程離散為[15，20] dρidt =jmj（vi-vj）·iWij

本文應用層流黏性對邊界處進行近似邊界層理論處理，更加接近實際物理情況，與傳統(tǒng)的單一排斥力等邊界處理方法相比，減緩了流體與壁面沖擊過程中的數(shù)值振蕩，最明顯的就是降低了沖擊時刻的壓力梯度，使得數(shù)值解更加光滑。該黏性模型既能保證線動量和角動量的守恒性，能量的耗散也較小，并且具有較好的層流模擬能力。在邊界層處選用該模型，并對邊界周圍的粒子運動進行修正，能最大程度降低黏性力模型造成的數(shù)值影響。另外，在流體計算過程中，邊界層內應用黏性方程（5）是為了更加接近邊界層理論，而對于主流區(qū)即非邊界層內的流體則采用了SPH常用的人工黏性來處理流體。

本文在進行邊界處理時，將接近壁面的流體視為層流，在耦合動力邊界附近引入層流黏性近似邊界層理論，有利于提高潰壩問題的數(shù)值模擬精確度，這將在下一節(jié)的模擬計算中得到驗證。

2?潰壩問題的數(shù)值仿真

為了與實驗更好地比較，本文計算的物理模型與文獻[2]的實驗模型相一致，如圖2所示。容器長度為1 610 mm，高度為600 mm，初始時水體的尺寸為600 mm×300 mm，水體寬度為600 mm，高度為300 mm。H1～H4是觀察水體實時表面高度的具體位置。

2.1?流體外部輪廓比較

為了更系統(tǒng)地和文獻[2]中的實驗結果進行比較，本文對流體外輪廓的高度和時間進行相同的無量綱化處理，無量綱液面高度的h*表達式為h*=h/H，無量綱時間t*的表達式為t=tg/H。此處h為實時液面高度，H為初始時刻的液面高度。圖3給出了不同時刻實驗與SPH計算的外部輪廓圖，以及不同時刻相對壓力的分布情況，其中時間后面括號內的數(shù)值是相對應時刻的無量綱時間。從圖中可以看出，在不同時刻，SPH計算結果的輪廓和實驗的輪廓基本吻合，尤其在t=5.851 6和t=6.671 1時，液體的飛濺現(xiàn)象明顯，也基本和實驗中飛濺現(xiàn)象的輪廓相一致。圖3中右側彩色云圖展示的是不同時刻流體的相對壓力的分布情況，從流場的壓力分布可以看出，在不同時刻，壓力分布光滑、連續(xù)。圖4所示的是t=4.931 0時刻右側拐角處的流場圖，從流場的具體分布可以看出，流場順滑，并且明顯地展現(xiàn)出流體流動方向改變而產生的渦，在邊界附近，流體的速度相對較小，表明了近似邊界層的修正起到了明顯作用。

結合圖3可以明顯的看出：在水流沖擊右側壁面后，隨著水流的翻卷并傾覆，與主流區(qū)域相融合之前，不同位置液面高度隨時間的變化與實驗吻合相當好。但當上涌的水流與主流區(qū)域相融合之后，即圖5中（b）～（d）發(fā)生明顯的拐點之后，SPH計算結果與實驗結果相對符合度差了些，但總體趨勢基本保持一致。主要是由于水-氣相互作用導致的能量耗散考慮不足，在自由表面附近，液體間的碰撞融合過程伴隨有黏性耗散作用，對自由表面的形成和自由表面的演化過程有重要影響，這是本文計算和實驗產生差異的重要原因。

2.3?監(jiān)測點壓力變化比較

圖6展示了右側壁面不同高度處壓力隨時間的變化過程，其中（a）～（d）分別展示了距容器底部3，15，30，80 mm高處的SPH計算壓力和100次傳感器測試壓力的95%置信區(qū)間的對比曲線，圖中顯示的壓力為無量綱壓力p，無量綱處理方式為實時壓力與最大靜水壓力的比值，即p=p/（ρ0gH）。在3 mm處，壓力幾乎全部落入置信區(qū)間內部，其中在流體剛沖擊到右側壁面時，p*出現(xiàn)了最大值，且p=6.79，大于置信區(qū)間的最大壓力值，這是因為流體與右側壁面發(fā)生撞擊，壓力從0在極短時間內快速增加，壓力梯度變化很大導致的。應用SPH方法研究文獻[2]中的潰壩水流模型成果很多[5，8-11]，本文與文獻[5]的結果做了對比，其中最大壓力與文獻[5]中的最大壓力值相差不大。在水流與右側壁面撞擊，經過翻卷并傾覆，與主流區(qū)域相融合后壓力也產生了波動，這主要是由于液體碰撞融合時，液面處于振蕩狀態(tài)，壓力波向周圍傳遞引起附近區(qū)域流體的壓力產生振蕩，這也是其余監(jiān)測點產生壓力波動的主要原因之一。

15，30 mm處的壓力大部分落入置信區(qū)間之內，相比于文獻[5]，大部分時間段的壓力值大于置信區(qū)間的最大值，本文的部分計算結果略低于置信區(qū)間的最小值，但都比較接近于置信區(qū)間的最小值，這主要是固壁附近流體粒子受邊界斥力和黏性力的共同作用后，固壁附近的流體粒子與壁面間的距離稍大于粒子間的初始間距，致使壓力值偏小。

在80 mm處，幾乎整個時間段域內的壓力值均小于置信區(qū)間的最小值，流體經過翻卷并傾覆，與主流區(qū)域相融合后，壓力開始升高，致使后半段時間內的壓力均落于實驗的置信區(qū)間之內，主要是因為加入了近似邊界層理論，流體與壁面間產生了相應的黏性效應，使得流體與壁面的接觸更加自然合理，從而使得后半部分的壓力升高，更接近實驗結果。

本文對不考慮層流邊界層影響的條件也進行了數(shù)值計算，發(fā)現(xiàn)不考慮邊界層效應時，壓力波動較大，特別是在流體沖擊壁面的瞬間，p*>10.0，在流體經過翻卷并傾覆后，流體的壓力波動幅度和振蕩時間都比考慮邊界層效應大很多，這也證明了文中方法是有效可行的。另外文中近似邊界層的厚度取為2.4mm，主要是從SPH邊界計算方法的角度進行考慮的，具體值的大小是根據(jù)邊界粒子的支持域半徑來定的，以后的工作可以結合邊界區(qū)域實時雷諾數(shù)的大小和自適應粒子分辨率技術，詳細研究邊界層厚度的選取對不同流動問題的影響。

3?結 論

本文在耦合動力邊界的基礎上，在壁面附近引入層流黏性來近似邊界層理論，應用改進的邊界處理方法對經典潰壩問題進行數(shù)值模擬，并與實驗結果進行了對比研究，得到以下結論。

（1） 較完整地模擬到了水流與壁面相互作用而產生的水花飛濺 、融合 、水流反彈、自由表面變化以及近壁面處水流的劇烈變形等多種復雜的物理現(xiàn)象，其外部輪廓與實驗非常吻合。

（2） 自由液面的高度在水流發(fā)生翻卷并傾覆與主流區(qū)域相融合之前相一致，在自由表面融合過程中，由于液面間的沖擊能量耗散，會導致融合后的液面高度存在一些差異。

（3） 數(shù)值模擬的壓力基本落在100次實驗研究的95%置信區(qū)間之內，近似邊界層理論在80 mm處的壓力變化的預測更加準確。

參考文獻：

[1]?李云，李君.潰壩模型試驗研究綜述[J].水科學進展，2009，20（2）：304-310.

[2]?Lobovsky L，Botia-Vera E，Castellana F，et al.Experimental investigation of dynamic pressure loads during dam break[J].Journal of Fluids and Structures，2014（48）：407-434.

[3]?Lee T，Zhou Z，Cao Y.Numerical simulations of hydraulic jumps in water sloshing and water impacting[J].Journal of Fluids Engineering，2002，124（1）：215-226.

[4]?Buchner B.Green water on ship-type offshore structures[D].Delft：Delft University of Technology，2002.

[5]?Cercos-Pita J L.AQUAgpusph，a new free 3D SPH solver accelerated with OpenCL[J].Computer Physics Communications，2015（192）：295-312.

[6]?Nguyen V T，Park W G.A free surface flow solver for complex three‐dimensional water impact problems based on the VOF method[J].International Journal for Numerical Methods in Fluids，2016，82（1）：3-34.

[7]?Mokrani C，Abadie S.Conditions for peak pressure stability in VOF simulations of dam break flow impact[J].Journal of Fluids and Structures，2016（62）：86-103.

[8]?Sun P， Ming F， Zhang A. Numerical simulation of interactions between free surface and rigid body using a robust SPH method[J].Ocean Engineering，2015（98）：32-49.

[9]?Aureli F，Dazzi S，Maranzoni A，et al.Experimental and numerical evaluation of the force due to the impact of a dam-break wave on a structure[J].Advances in Water Resources，2015（76）：29-42.

[10]?Sun P，Ming F，Zhang A.Numerical simulation of interactions between free surface and rigid body using a robust SPH method[J].Ocean Engineering，2015（98）：32-49.

[11]?Oger G，Marrone S，Le Touzé D，et al.SPH accuracy improvement through the combination of a quasi-Lagrangian shifting transport velocity and consistent ALE formalisms[J].Journal of Computational Physics，2016（313）：76-98.

[12]?彭榮強.應用等位函數(shù)法模擬潰壩流動[J].長江科學院院報，2008，25（2）：11-12.

[13]?Violeau D，Rogers B D.Smoothed particle hydrodynamics （SPH） for free-surface flows： past， present and future[J].Journal of Hydraulic Research，2016，54（1）：1-26.

[14]?Shadloo M S，Oger G，Le Touzé D.Smoothed particle hydrodynamics method for fluid flows，towards industrial applications：Motivations， current state， and challenges[J].Computers & Fluids，2016（136）：11-34.

[15]?Liu G R，Liu M B.Smoothed particle hydrodynamics： a meshfree particle method[M].Newyork：World Scientific，2003.

[16]?Liu M B，Liu G R.Smoothed particle hydrodynamics （SPH）：an overview and recent developments[J].Archives of Computational Methods in Engineering，2010，17（1）：25-76.

[17]?Yue W，Lin C L，Patel V C.Numerical simulation of unsteady multidimensional free surface motions by level set method[J].International Journal for Numerical Methods in Fluids，2003，42（8）：853-884.

[18]?Lv X，Zou Q，Zhao Y，et al.A novel coupled level set and volume of fluid method for sharp interface capturing on 3D tetrahedral grids[J].Journal of Computational Physics，2010，229（7）：2573-2604.

[19]?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（3）：375-389.

[20]?Monaghan J J.Smoothed particle hydrodynamics[J].Annual Review of Astronomy and Astrophysics，1992，30（1）：543-574.

[21]?Shao J R，Li H Q，Liu G R，et al.An improved SPH method for modeling liquid sloshing dynamics[J].Computers & Structures，2012（100）：18-26.

[22]?Macià F，Sánchez J M，Souto‐Iglesias A，et al.WCSPH viscosity diffusion processes in vortex flows[J].International Journal for Numerical Methods in Fluids，2012，69（3）：509-533.

[23]?Colagrossi A，Antuono M，Le Touzé D.Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model[J].Physical Review E，2009，79（5）：056701.

[24]?Monaghan J J，Gingold R A.Shock simulation by the particle method SPH[J].Journal of Computational Physics，1983，52（2）：374-389.

（編輯：胡旭東）