邵 晨，黃劍峰

（云南農(nóng)業(yè)大學(xué)水利學(xué)院，昆明650201）

0 引言

潰壩是一種災(zāi)害性的水流現(xiàn)象，其洪峰流量、運動速度、破壞力遠(yuǎn)大于一般暴雨洪水或融雪洪水。除了其超標(biāo)準(zhǔn)洪水、冰凌以及地震等致大壩潰壩的自然因素外，還有設(shè)計不周、施工不良以及管理不善等人為因素。潰壩的破壞力會帶來巨大損失［1］，所以潰壩數(shù)值模擬就顯得尤為重要，模擬結(jié)果不僅可以分析水流的流動過程，還可以用于推演壩前設(shè)置不同障礙物的消能效果，而且能夠?qū)紊鼡p失進(jìn)行估算［2］。

計算流體力學(xué)（CFD）可以在虛擬環(huán)境中仿真分析流體流動情況，降低了時間和成本，是一種可以代替實驗研究的手段。傳統(tǒng)的求解Navier-Stokes 方程的CFD 解算器一般都是基于網(wǎng)格的求解方法，因此需要經(jīng)過一個耗時的網(wǎng)格劃分過程，網(wǎng)格質(zhì)量不僅制約著數(shù)值解的穩(wěn)定性，而且影響著解的質(zhì)量，且若干等尺寸的網(wǎng)格會提供不同的解，這會給解帶來不確定性。大多數(shù)傳統(tǒng)的編碼通過加入人為的數(shù)值黏度來確定解，這種黏度會很大程度上改變和抑制解的穩(wěn)定性，并且對非恒定流和聲學(xué)的精度有很大的負(fù)面影響。如果不考慮一些簡化和假設(shè)，那么Navier-Stokes方程是沒有解析解的，而且由于計算的限制，它們的直接數(shù)值模擬（Direct numerical simulation）在實際工程中幾乎是不可能的，所以僅限于學(xué)術(shù)應(yīng)用，因此需要建立模型。然而，像雷諾平均方程（RANS）湍流模型都是經(jīng)驗?zāi)Ｐ?，它們的缺點是依賴大量常數(shù)，為了正確地模擬湍流就要必須仔細(xì)校準(zhǔn)這些常數(shù)。發(fā)展的大渦模擬（LES）湍流模型因其與使用Navier-Stokes 求解器的RANS 模型相比計算成本較高，所以并沒有得到廣泛的應(yīng)用。LBM 是一種介于宏觀與微觀之間的介觀方法，它不考慮單個粒子的運動，而是將所有粒子的運動視為一個整體，它同時具有微觀方法和宏觀方法的優(yōu)勢，其更擅長處理復(fù)雜邊界，具有良好的并行性，且其獲取壓力場的方式不像傳統(tǒng)N-S 方程需要大量計算求解泊松方程，LBM 是通過計算流體狀態(tài)方程獲得壓力場。通過LBM 模型在處理破碎交界面的優(yōu)勢與VOF相結(jié)合，可以更好的觀察到液體的飛濺及破碎。

為了對潰壩后水流的流動過程進(jìn)行數(shù)值計算，開發(fā)出如SPH 法［3］等二維模型，該方法能簡便、快速、有效地進(jìn)行數(shù)值模擬。寧聰?shù)龋?］利用二維模型對不同閘門開度下潰壩水流的淹沒范圍及最大流速進(jìn)行了分析；劉慧玲等［5］利用SPH 方法對潰壩水流進(jìn)行模擬，通過實驗與數(shù)值模擬的潰壩水流外部輪廓進(jìn)行比較，驗證了SPH方法的準(zhǔn)確性；SWEGLE等［6］發(fā)現(xiàn)SPH法存在著如自由表面的計算精度不夠準(zhǔn)確，計算應(yīng)力時會出現(xiàn)計算不穩(wěn)定等現(xiàn)象；陳薇等［7］根據(jù)改進(jìn)的核函數(shù)和改進(jìn)的動量方程的離散形式相結(jié)合后驗證發(fā)現(xiàn)SPH 法穩(wěn)定性和模擬效率都有明顯提高；繆吉倫等［8］論述了SPH 法的基本理論并利用該方法模擬了潰壩水流沖擊障礙物的飛濺及融合。對于進(jìn)一步開發(fā)的三維模型，林長強(qiáng)［9］利用FLUENT 對土石壩逐漸潰壩水流模擬；曾丹等［10］在基于Flow-3D，采用RNGk-ε模型和VOF 法，建立了分析潰壩洪水的水流三維數(shù)學(xué)模型，模擬彎道潰壩洪水的演進(jìn)過程，其計算精度與水流模擬的流動過程相對較好，但其又存在著需要復(fù)雜且耗時的網(wǎng)格劃分等問題。近年來，基于最小動力學(xué)模型的Boltzmann 方程越來越受歡迎。LBM 方法最初是用于去除統(tǒng)計噪聲并實現(xiàn)更好的伽利略不變性［11，12］，由于其與動力學(xué)理論的密切聯(lián)系所提供的靈活性，LBM 方法可用于模擬多種物理現(xiàn)象，張建民等［13］對LBM 應(yīng)用在多相流中的研究現(xiàn)狀進(jìn)行了綜述，旨在進(jìn)一步擴(kuò)大LBM 方法在多相流研究中的適用性。本文將采用LBM 方法對瞬時全潰壩和局部瞬時潰壩進(jìn)行數(shù)值模擬，并與前人的模擬實驗結(jié)果進(jìn)行比較，利用LBM 方法對瞬時潰壩時在壩前設(shè)置不同消能坎對其消能能力進(jìn)行比較，并分析水流對障礙物的作用力，從而表明基于LBM 方法在潰壩數(shù)值模擬中具有很大的使用價值。

1 數(shù)值模擬方法

1.1 格子玻爾茲曼方法

格子Boltzmann 方法（LBM）是從格子氣自動機(jī)（LGA）方法發(fā)展來的，與LGA 方案所使用的布爾邏輯不同的是，LBM 利用具有實變量的統(tǒng)計分布函數(shù)fi，通過構(gòu)建質(zhì)量守恒和線性動量來計算數(shù)據(jù)。

定義格子Boltzmann的輸運方程如下：

式中：fi是方向i上的粒子分布函數(shù)；ei是相應(yīng)的離散速度；Ωi是碰撞算子。

LBM 的撞擊流方法可解釋為連續(xù)Boltzmann 方程的離散化。撞擊流或傳播步驟模擬的是粒子分布函數(shù)沿離散方向的平流，而大多數(shù)模型是由碰撞算子模擬的，碰撞算子對模擬的數(shù)值穩(wěn)定性也有很大影響。最常見的方法是使用基于Bhatna?gar-Gross-Krook（BGK）近似的單弛豫時間（SRT），其定義如下所示：

式中：τ是松弛時間參數(shù)，與宏觀黏度v相關(guān)［14］，如下所示：

這里的cs是格子聲速，p為流體密度，u是宏觀速度，δ是克羅內(nèi)克符號，α和β表示方程中出現(xiàn)的向量的不同空間分量，wi是為了保持各向同性而建立的加權(quán)常數(shù)，對于常用的三維格子離散速度模型D3Q19模型（見圖1），有：

圖1 D3Q19格子模型Fig.1 D3Q19 lattice model

通過Chapman-Enskog 展開，可以證明所得到的結(jié)果可以重現(xiàn)低馬赫數(shù)下的水動力狀態(tài)［15，16］。單弛豫時間法因其簡單而常被使用，但它不適用于高馬赫數(shù)的應(yīng)用，而且容易出現(xiàn)數(shù)值不穩(wěn)定。用多弛豫時間（MRT）碰撞算子解決了BGK 的一些限制，其碰撞過程在矩形空間中而不是在通常所用的速度空間中進(jìn)行，其定于如下：

另一種旨在克服BGK 方法局限性的方法是熵晶格Boltzmann（ELBM）方法，它可以依賴于單個弛豫時間，然而從計算的角度來看，這種方法特別昂貴，因此沒有在實際工程應(yīng)用中使用。中心力矩格子Boltzmann 方法中的碰撞算子是基于多重弛豫時間的方法。然而，與標(biāo)準(zhǔn)MRT不同的是散射算子是在中心矩空間中實現(xiàn)的。松弛過程是在一個移動的參考系中通過改變離散粒子的速度和局部宏觀速度來實現(xiàn)的，這自然地提高了給定速度集的伽利略不變性和數(shù)值穩(wěn)定性［19］。原始力矩可以定義為：

中心力矩定義為：

式中：ux，uy和uz是宏觀速度的組成部分；μ為原始力矩；k、l和m分別為x、y和z方向上的力矩階數(shù)。

1.2 湍流模型

湍流建模采用的方法是大渦模擬（LES）。該方法引入了一種附加的黏度系數(shù)ut，稱為紊流渦流黏度Vt，以模擬亞網(wǎng)格湍流。本文使用的LES 模型是壁面適應(yīng)局部渦流（WALE）黏度模型，該模型提供了一致的局部渦流黏度和近壁面行為。具體定義如下：

式中：?f=Cw?x是濾波尺度；?x為單位網(wǎng)格尺度；S是分解尺度的應(yīng)變率張量；常數(shù)Cw通常為0.325；u為動力黏度系數(shù)；δαβ為克羅內(nèi)克爾符號；α，β，r為張量下標(biāo)。

式中：y為與壁面的法向距離；x為與壁面相切的局部流動方向；uτ為表面摩擦速度；τw為湍流壁面剪應(yīng)力；dpw/dx為壁面壓力梯度的特征速度；u為離壁面給定距離處的平均速度；f1和f2為插值函數(shù)。

邊界層的速度場是通過y+得到的；y+取決于第一個晶格從幾何體壁y到第一個晶格的速度uc之間的距離。

1.3 VOF法

VOF 法是根據(jù)不同時刻單元網(wǎng)格內(nèi)流體容積的改變實現(xiàn)對自由表面的追蹤，其具體定義可描述為：假設(shè)存在某一流體ρ，某一時刻時流體ρ所占單元網(wǎng)格的體積函數(shù)為ɑ，當(dāng)ɑ=1 時，則單元網(wǎng)格內(nèi)充滿流體ρ；當(dāng)ɑ=0 時，該單元網(wǎng)格內(nèi)無流體ρ；當(dāng)0<ɑ<1時，則說明該單元處于氣體與流體ρ交界面。

2 模型驗證

本文采用所建立的三維模型對全潰壩和局部潰壩2個物理模型進(jìn)行數(shù)值模擬，并將計算結(jié)果與楊哲豪等［21］數(shù)值模擬結(jié)果進(jìn)行比較，驗證分析模型的準(zhǔn)確性。本文中默認(rèn)取水的密度為ρ=1 000 kg/m3，采用牛頓黏度模型，動力黏度為0.001 Pa?s，重力加速度為9.81 m/s2。

2.1 瞬時全潰壩

參考Al-Faesly等［22］幾何模型，建立封閉的幾何箱體進(jìn)行數(shù)值模擬。瞬時全潰壩數(shù)值模擬俯視圖如圖2所示，采用長5.6 m，寬2.7 m，高為1 m 的矩形箱體作為水壩儲水，下游水槽設(shè)置為長為14 m，寬為1.32 m，高為1 m 的矩形箱體，因?qū)Ρ葘嶒灢捎玫氖窍掠文┒怂蹫殚_邊界，為了保證實驗的準(zhǔn)確性對下游水槽的封閉邊界進(jìn)行了適當(dāng)加長。此外，實驗還設(shè)置了兩個傳感器來監(jiān)測潰壩時洪水的高度。水壩的初始水位設(shè)為0.85 m，模擬時間為8 s，模擬精度為0.05 m，潰壩方式為瞬時潰壩。

圖2 全局瞬時潰壩數(shù)值模擬俯視圖（單位：m）Fig.2 Top view of numerical simulation of global instantaneous dam break

通過數(shù)值模擬水流的流動過程，結(jié)果如圖3所示。水體在初始位置時開始向下游水槽流動，在重力加速度的作用下，潰壩水流在流動過程中的速度形成了非常鮮明的分層現(xiàn)象。從圖3（a）看到水流到達(dá)了傳感器1 的位置，從圖3（b）看到在3.16 s洪水到達(dá)左壁面，由圖3（c）看到在4.12 s時擊打左壁面后水流飛濺形成一定高度的水舌。

圖3 瞬時全潰壩水流的流動過程Fig.3 The flow process of global instantaneous dam break

在傳感器1 監(jiān)測潰壩時洪水的水位高度變化，通過記錄水位高度變化的并與文獻(xiàn)［21］進(jìn)行比較，對比結(jié)果如圖4（a）所示。在本文所建立的三維模型下進(jìn)行數(shù)值模擬，計算結(jié)果到達(dá)傳感器1監(jiān)測點位置的時間與文獻(xiàn)［21］實驗到達(dá)監(jiān)測點的時間更相近，二維數(shù)值模擬下水流到達(dá)傳感器1 的時間略早于實驗到達(dá)的時間，由于實驗中水泵一直運作導(dǎo)致模擬過程中水位下降比實驗更快，所以模擬結(jié)果存在一定的偏差。圖4（b）中8 s內(nèi)傳感器2監(jiān)測點的水位變化計算結(jié)果與文獻(xiàn)［21］計算結(jié)果基本相同，可以看出本文的模擬結(jié)果更接近實驗中水位下降高度的監(jiān)測數(shù)據(jù)。通過圖4 數(shù)據(jù)結(jié)果表明基于LBM 的方法在潰壩數(shù)值模擬計算中是準(zhǔn)確的。

圖4 相同時間內(nèi)水位變化比較圖Fig.4 Comparison of water level changes in the same time

2.2 局部瞬時潰壩

上文驗證了瞬時全潰壩的準(zhǔn)確性，為了進(jìn)一步驗證模型在局部瞬時潰壩時的準(zhǔn)確性，本文根據(jù)Aureli 等［23］的局部潰壩物理模型，建立封閉的幾何箱體進(jìn)行數(shù)值模擬。局部潰壩的數(shù)值模擬俯視圖如圖5所示，設(shè)立一個長2.6 m，寬1.2 m，高1 m 的封閉箱體，設(shè)水庫初始水位為0.1 m，水庫長0.8 m，寬1.2 m，在水壩中央位置設(shè)置0.3 m 寬大壩以達(dá)到局部瞬時潰壩的目的。在離水壩0.5 m 下游處設(shè)置一長0.155 m，寬0.3 m，高0.2 m 的矩形箱體障礙物，并在適當(dāng)位置設(shè)置兩個速度傳感器，模擬時間為1.5 s，模擬精度為0.009 m。

圖5 局部潰壩的數(shù)值模擬俯視圖（單位：m）Fig.5 Top view of numerical simulation of local dam break

通過LBM 方法對局部潰壩進(jìn)行數(shù)值模擬并選取四個不同時刻，與文獻(xiàn)［21］現(xiàn)場實驗和計算水深云圖進(jìn)行比較，如圖6所示。圖6（b）可以看出在0.38 s 水流到達(dá)障礙物，圖6（c）看出在水流接觸到障礙物位置的水深明顯爬升，在圖6（d）中水流繞障礙物向兩側(cè)流動。在LBM 方法的數(shù)值模擬中，VOF下呈現(xiàn)的模擬結(jié)果與文獻(xiàn)［21］實驗拍攝圖形與計算云圖吻合較好，表明LBM方法可以準(zhǔn)確的模擬局部瞬時潰壩水流的流動情況。

黨的十八大報告提出，要實施創(chuàng)新驅(qū)動發(fā)展戰(zhàn)略，強(qiáng)調(diào)科技創(chuàng)新是提高社會生產(chǎn)力和綜合國力的戰(zhàn)略支撐，必須擺在國家發(fā)展全局的核心位置。這充分表明了我們黨依靠創(chuàng)新實現(xiàn)經(jīng)濟(jì)持續(xù)健康發(fā)展的堅定決心和對科技創(chuàng)新的高度重視。當(dāng)前，東營市正處于全面實施黃藍(lán)國家戰(zhàn)略，加快推進(jìn)“率先全面建成小康社會、率先建成生態(tài)文明典范城市”（兩個率先）建設(shè)目標(biāo)的關(guān)鍵時期，對科技進(jìn)步和創(chuàng)新提出了更加全面、更加緊迫的需求，必須深刻認(rèn)識并充分發(fā)揮科技創(chuàng)新的支撐引領(lǐng)作用，制定和實施創(chuàng)新驅(qū)動戰(zhàn)略，完善區(qū)域科技創(chuàng)新體系，全力依靠創(chuàng)新驅(qū)動提高發(fā)展的質(zhì)量和效益。

圖6 局部潰壩數(shù)值模擬對比圖Fig.6 Comparison of numerical simulation of local dam break

為了比較不同位置的速度變化，對不同位置傳感器的速度進(jìn)行計算，結(jié)果如圖7所示。由圖7（a）傳感器速度的變化可以看出，潰壩水流在0.38 s左右到達(dá)傳感器位置，流速陡然上升到0.155 m/s，隨后水流到達(dá)障礙物后速度急劇下降，這是由于障礙物的阻礙作用；在0.6～0.78 s 流速隨著潰壩水流演進(jìn)流速加快；在0.78～0.96 s 隨著障礙物上游水深爬升，流速逐漸變慢；0.96～1.04 s 障礙物上游儲存的水開始向兩側(cè)流動，流速隨之增大。圖7（b）可以看出水流在0.52 s 到達(dá)傳感器4 流速陡然上升到0.85 m/s，在1.08 s流速達(dá)到最大；當(dāng)潰壩水流撞擊障礙物后，水流從障礙物的兩側(cè)向下游流動，通過比較可以看出圖7（b）水流的速度大于圖7（a），這和實際中結(jié)構(gòu)兩側(cè)水流速度快相符。

圖7 不同位置傳感器速度變化圖Fig.7 Velocity variation of sensors at different positions

為了進(jìn)一步研究局部潰壩水流對障礙物的作用力，對1.5 s內(nèi)水流對障礙物的作用力進(jìn)行計算，結(jié)果如圖8所示。水流在0.38～0.9 s撞擊障礙物然后繞障礙物開始流動，因水流對障礙物的不斷撞擊，其作用力變化劇烈；水流在0.9～1.5 s 已包裹住整個障礙物，水流的流動也趨于平緩，水流對障礙物的作用力變化劇烈程度變小。水流對障礙物沖刷的作用力變化劇烈程度體現(xiàn)了水流流動的不穩(wěn)定性，而在文獻(xiàn)［21］中利用二維數(shù)值模擬所計算出的潰壩洪水對障礙物的沖擊荷載曲線圖相對平滑，不能體現(xiàn)出由于洪水對障礙物的劇烈沖刷所帶來的沖擊荷載的劇烈變化程度。

圖8 水流對障礙物的作用力Fig.8 Force of flow on Obstacles

3 兩種消能坎的潰壩流動模擬

通過與前人實驗比較后可證明基于格子Boltzmann 方法的準(zhǔn)確性，利用該方法分別模擬無孔障礙物和有孔障礙物兩種情形下的自由表面流動現(xiàn)象。參考劉維平等［24］幾何模型，在離封閉箱體左側(cè)0.6 m 處布設(shè)了一個障礙物，通過比較兩種模型中水流撞擊障礙物及壁面時的壓強(qiáng)，分析兩種障礙物的消能效果。

無孔障礙物與有孔障礙物模型側(cè)視圖模型布置如圖9所示，建立一個長1.2 m，寬0.8 m，高1.2 m 的封閉箱體，在箱體左壁離右側(cè)0.6 m 處設(shè)置一個長0.1 m，寬0.8 m，高0.2 m 的矩形無孔障礙物，如圖9（a）所示。在其他條件相同的情況下，在障礙物上設(shè)長0.1 m，寬0.8 m，高0.06 m 的孔徑，圖9（b）所示。水體密度為ρ=1 000 kg/m3，采用牛頓黏度模型，取用動力黏度為0.001 Pa?s，重力加速度為9.81 m/s2，潰壩方式為瞬時潰壩，模擬時間為0.9 s，模擬精度為0.01 m。

圖9 無孔障礙物與有孔障礙物模型側(cè)視圖（單位：m）Fig.9 Side view of nonporous and perforated obstacle models

相同時刻兩種不同消能方式的水流流動運動情況如圖10所示。在瞬時潰壩后水體發(fā)生移動，水體撞擊障礙物后向上飛濺。在t=0.14 s時，水流到達(dá)障礙物，在0.3 s時明顯可觀察到水流拍打障礙物激起水浪。無孔障礙物水舌飛濺高度最高達(dá)到1.2 m，有孔障礙物的水舌飛濺高度達(dá)1 m，速度的變化會引起數(shù)值模擬中流體顏色的不同，可以觀察到無孔障礙物水流速度更快，加之飛濺后水舌的高度過高，考慮到下落過程重力勢能的影響會導(dǎo)致其破壞性更強(qiáng)。

圖10 相同時刻兩種消能方式水流運動對比Fig.10 Comparison of flow motion of two energy dissipation methods at the same time

選取2個典型時刻兩種消能坎的VOF 等值面如圖11所示。從圖11（a）可以看到水流在撞擊障礙物飛濺到一定高度后波浪發(fā)生破碎，從圖11（b）和圖11（c）看到在水流撞擊壁面后波浪再一次發(fā)生破碎，通過等值面中水流躍起形成的水舌高度以及波浪的破碎程度看出有孔障礙物消能效果更明顯。

圖11 兩種消能坎的VOF等值面圖Fig.11 VOF isosurface of two kinds of energy dissipation sills

為了驗證經(jīng)過有孔障礙物與無孔障礙物時兩種消能坎的消能效果，分別在壁面不同高度設(shè)置了傳感器監(jiān)測水流沖擊壁面的壓強(qiáng)，水體撞擊壁面的壓強(qiáng)比較見圖12所示。通過比較可見無孔障礙物產(chǎn)生的壓強(qiáng)更大，因此可在消能障礙物上設(shè)置一定大小的孔隙，保證障礙物有更好的消能效果。

圖12 水體撞擊左壁面的壓強(qiáng)比較Fig.12 Pressure comparison of water impact on left wall

為了比較水流對障礙物的作用力，對水流沖擊障礙物時的作用力進(jìn)行計算，結(jié)果如圖13所示。分別對水流在有孔障礙物的上部結(jié)構(gòu)與下部結(jié)構(gòu)及無孔障礙物的作用力進(jìn)行數(shù)值模擬計算，水流在0.14 s到達(dá)障礙物時，有孔障礙物的下部結(jié)構(gòu)與無孔障礙物受到水流的作用力基本相同，無孔障礙物在0.2 s左右受到水流作用力達(dá)到最大，最大值為1 589 N。高于障礙物的水流在潰壩時由于重力加速度的作用，水流對障礙物作用力更大，在高于障礙物水流沖刷過后，水流對障礙物作用力的波動頻率趨于穩(wěn)定。根據(jù)圖13的作用力對比結(jié)果，在相同的潰壩情況下，有孔障礙物抵抗洪水沖刷所需要的強(qiáng)度要求遠(yuǎn)小于無孔障礙物，無孔障礙物對結(jié)構(gòu)的強(qiáng)度要求更高。

圖13 水流對障礙物的作用力Fig.13 Force of flow on Obstacles

4 結(jié)論

本文在格子Boltzmann方法的基礎(chǔ)上，采用LES模型和VOF法，建立了分析了潰壩水流分析的三維模型，通過一系列數(shù)值模擬得到如下結(jié)論。

（1）驗證了該模型在瞬時全潰壩水流流動過程的準(zhǔn)確性，本文方法與實驗結(jié)果更相近在，在三維模擬下可以看到自由面流動情況，4.12 s 到達(dá)封閉箱體撞擊左壁后可看到濺起一定高度的水舌。

（2）驗證了該模型在局部瞬時潰壩水流流動過程的準(zhǔn)確性，可以看到水流到達(dá)障礙物時，由于受到障礙物的阻擋，障礙物上游水流高度向上爬升，至一定高度后回落。隨后，水流從障礙物兩側(cè)向下游運動；水流在向兩側(cè)變窄的河道流動流速加快，隨著水流繞障礙物兩側(cè)向下游流動，水流對障礙物的作用力波動頻率變小。

（3）通過模擬了兩種不同消能坎的潰壩水流的流動情況，可以看到潰壩洪水沖擊障礙物后水流飛濺的高度以及流動情況，通過壓強(qiáng)的比較得出了有孔障礙物做消能坎時消能效果更好；通過VOF等值面的水流流動情況可以看到水流撞擊無孔障礙物后波浪破碎程度更明顯；有孔障礙物比無孔障礙物受到水流作用力削減了60%，水流對無孔障礙物沖刷的作用力比有孔障礙物更大。□