劉慧玲 李海橋

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

關(guān)?鍵?詞：潰壩水流模擬; 自由液面; 耦合動(dòng)力邊界; 光滑粒子流體動(dòng)力學(xué)（SPH）

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

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

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

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

1?光滑粒子流體動(dòng)力學(xué)

1.1?基本原理

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

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

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

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

為了與實(shí)驗(yàn)更好地比較，本文計(jì)算的物理模型與文獻(xiàn)[2]的實(shí)驗(yàn)?zāi)Ｐ拖嘁恢?，如圖2所示。容器長(zhǎng)度為1 610 mm，高度為600 mm，初始時(shí)水體的尺寸為600 mm×300 mm，水體寬度為600 mm，高度為300 mm。H1～H4是觀察水體實(shí)時(shí)表面高度的具體位置。

2.1?流體外部輪廓比較

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

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

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

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

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

在80 mm處，幾乎整個(gè)時(shí)間段域內(nèi)的壓力值均小于置信區(qū)間的最小值，流體經(jīng)過(guò)翻卷并傾覆，與主流區(qū)域相融合后，壓力開始升高，致使后半段時(shí)間內(nèi)的壓力均落于實(shí)驗(yàn)的置信區(qū)間之內(nèi)，主要是因?yàn)榧尤肓私七吔鐚永碚摚黧w與壁面間產(chǎn)生了相應(yīng)的黏性效應(yīng)，使得流體與壁面的接觸更加自然合理，從而使得后半部分的壓力升高，更接近實(shí)驗(yàn)結(jié)果。

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

3?結(jié) 論

本文在耦合動(dòng)力邊界的基礎(chǔ)上，在壁面附近引入層流黏性來(lái)近似邊界層理論，應(yīng)用改進(jìn)的邊界處理方法對(duì)經(jīng)典潰壩問(wèn)題進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比研究，得到以下結(jié)論。

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

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

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

參考文獻(xiàn)：

[1]?李云，李君.潰壩模型試驗(yàn)研究綜述[J].水科學(xué)進(jìn)展，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]?彭榮強(qiáng).應(yīng)用等位函數(shù)法模擬潰壩流動(dòng)[J].長(zhǎng)江科學(xué)院院報(bào)，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.

（編輯：胡旭東）