基于VOF 方法和DFBI 模型的潰壩水流數(shù)值模擬

馬云飛，李濤，熊進標

上海交通大學 機械與動力工程學院，上海 200240

潰壩問題在計算流體力學(CFD)中由于其自由表面變化多樣，常被用來驗證數(shù)值模擬的可靠性。流體體積函數(shù)(VOF)方法和動態(tài)流體-剛體相互作用(DFBI)模型的結合多用于船舶工業(yè)設計等領域，其在模擬液體自由表面對船舶作用產生的多自由度運動問題中表現(xiàn)良好，基于此本文討論了VOF 方法結合DFBI 模型處理潰壩沖擊問題的數(shù)值模擬方法。

近些年來，國內對于潰壩問題的數(shù)值模擬方法進行了大量研究。其中，張力方等[1]將光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)和格子玻爾茲曼方法(lattice Boltzmann method，LBM)應用于潰壩問題并進行對比研究，結果表明：SPH方法在自由面捕捉方面有優(yōu)勢，LBM 在水位、壓力等參數(shù)方面預測效果較好。王興華等[2]針對潰壩過程中水位變化和潰壩負波對潰壩進程的影響進行研究，使用VOF 方法對潰壩水面流動狀態(tài)進行數(shù)值模擬并與已有數(shù)據(jù)進行了比對驗證。馬利平等[3]采用源項法耦合了潰口演變模型DB-IWHR與基于GPU 加速技術的二維水動力模型，他們建立的耦合模型所得結果與實測吻合較好。陳本毅等[4]基于自主研發(fā)的緊致插值曲線數(shù)學模型研究了潰壩水槽閘門不同抽離方向的影響。曹引等[5]基于結構化網(wǎng)格建立了有限體積模型Hydro M2D，其在復雜地形和不規(guī)則邊界情境下可以較好地模擬潰壩問題，并就此進行了各方面驗證。王笑等[6]、王福良等[7]、董俊哲等[8]都對潰壩數(shù)值模擬方向有一定研究成果?？梢妼τ跐螁栴}的數(shù)值模擬研究成果十分豐碩，但是關于潰壩沖擊效應問題的數(shù)值模擬方法卻較為匱乏。

張力方等[1]研究了SPH 方法和LBM 方法對可運動障礙物位移的模擬，2 種方法均有一定效果。楊哲豪等[9]采用有限體積法和中心迎風格式建立了二位數(shù)值模型，該模型在模擬全局潰壩工況方面有良好表現(xiàn)，但卻無法模擬實際水流紊動過程。

基于此，本文使用商用CFD 軟件Star-CCM+的Overset 網(wǎng)格技術，結合VOF 和DFBI 模型模擬了擋板潰壩問題，還設計了擋板潰壩實驗，運用高速可視化裝置測量了潰壩液位變化及擋板擺角變化，為進一步驗證該數(shù)值模擬方法，結合實驗數(shù)據(jù)進行了研究分析。

1 控制方程

1.1 基本控制方程

擋板開啟過程滿足Navier-Stokes 方程，即

式中：v為t時刻的速度矢量，ρ為密度，p為壓力，γ為運動黏度。

1.2 DFBI 模型

為了準確模擬擋板的運動，需要運用DFBI 方法。DFBI 模型用于模擬剛體響應流體施加的壓力和剪切力，其控制方程為

式中：T為擋板的凈扭矩；ω是關于旋轉軸的角速度；I為關于旋轉軸的慣性矩，其展開式為

扭矩T主要取決于水流的壓力扭矩Tp和剪切力轉矩Tτ，滿足

式中：f為水流所作用的表面，pf為作用在表面f上的壓力，τf為作用在表面f上的剪切力，rf為擋板質心到面f中心的距離向量，αf為面f的面積向量。

考慮到實驗裝置尺寸不小，僅擋板就有近1 kg，摩擦力矩可以忽略不計，認為摩擦轉矩為0。DFBI 模型僅以建立剛體質心的運動來得到整個剛體的運動，簡單、有效且方便，但DFBI 模型不適用于剛體質心保持不變的情況。

1.3 VOF 方法

VOF 方法[10]是通過研究網(wǎng)格單元中流體和網(wǎng)格體積比函數(shù)來確定自由面，追蹤流體的變化，而不是追蹤自由液面上質點的運動。其在網(wǎng)格上該相流體質點存在的定義函數(shù)如式(2)，質量守恒形式如式(3)，差分形式如式(4)，差分網(wǎng)格如圖1 所示。

圖1 網(wǎng)格單元上的變量表示[6]

VOF 方法形式簡單，結果有效，在經(jīng)過許多發(fā)展和改進后能處理越來越復雜的問題，自由面模擬結果也越來越精確，在多個不同領域發(fā)揮著重要作用。

2 實驗設計

擋板潰壩實驗裝置的幾何尺寸如圖2 所示。實驗裝置通體為有機玻璃所制，左側水箱內尺寸為300 mm×300 mm×300 mm，右側蓄水池內尺寸為290 mm×350 mm×600 mm，擋板尺寸為320 mm×320 mm×10 mm，質量為1.21 kg。實驗中，采用高速相機(Phantom V710)在擋板側面和擋板正面分別進行觀測。拍攝過程頻率為500 Hz，圖片分辨率為1 280 px×800 px。在水箱左側外表面設置了標尺以確定最高液位。對于擋板擺角，取擋板外表面末端像素點作為參考點，根據(jù)其在運動過程中關于初始狀態(tài)擋板外表面和經(jīng)過外表面軸端點水平線的相對垂直距離得到擺角兩直角邊，使用反正切三角函數(shù)得到擺角值。本文就初始液位為100、150、200 mm 的3 個工況開展了觀測實驗。

圖2 實驗裝置

3 網(wǎng)格與物理模型

整個計算域由Overset 網(wǎng)格區(qū)域和背景網(wǎng)格區(qū)域組成：左側背景網(wǎng)格區(qū)域尺寸為300 mm×300 mm×300 mm；右側為Overset 網(wǎng)格運動的背景網(wǎng)格區(qū)域，其面積大于Overset 網(wǎng)格運動過程中掃過的區(qū)域面積，尺寸為290 mm×350 mm×325 mm；擋板垂直置于左側箱體側面，且高于箱體邊緣10 mm。

左側背景網(wǎng)格區(qū)域上表面和右側背景網(wǎng)格區(qū)域下表面設為壓力出口邊界。在擋板周邊設置Overset 邊界，形成外尺寸為350 mm×340 mm×50 mm的Overset 網(wǎng)格區(qū)域，并與背景區(qū)域之間形成Overset 網(wǎng)格界面，使用緊密接近和代替孔切削設置。圖3 為以Y軸為法向經(jīng)過擋板幾何中點截面的網(wǎng)格場景。Overset 網(wǎng)格與背景網(wǎng)格均以局部結構化網(wǎng)格處理，采用切割體網(wǎng)格外加棱柱層，網(wǎng)格尺寸為7 mm，注意Overset 網(wǎng)格尺寸需與背景網(wǎng)格尺寸盡量保持一致，否則在受體網(wǎng)格單元層的數(shù)據(jù)插值精度將受到影響。對潰壩流體主要活動區(qū)和擋板運動掠掃區(qū)域這2 個重點區(qū)域特意進行了局部網(wǎng)格加密。使用體積控制加密，其網(wǎng)格絕對尺寸為5 mm 與2.5 mm。生成背景區(qū)域85.8 萬網(wǎng)格，重疊區(qū)域16.2 萬網(wǎng)格。

圖3 網(wǎng)格與幾何條件

物理模型采用可實現(xiàn)的雙層k-ε湍流模型與VOF 流體域體積方法，外加重力影響。擋板設置DFBI 單自由度旋轉運動，旋轉軸設于擋板上表面中線，方向為-Y軸方向，質心為擋板幾何中心，慣性矩對角分量根據(jù)式(1)計算。

隱式不定常設置為一階時間離散，由于受體網(wǎng)格單元層會根據(jù)剛體運動更新邊界，原則上在一個時間步長內剛體運動最大距離為網(wǎng)格單元尺寸。設置時間步長為0.002 s。為了每個時間步長內部迭代收斂完全，將最大內部迭代次數(shù)設置為20，最大物理時間設置為3 s。

4 模擬結果與分析

4.1 網(wǎng)格敏感性

為了分析網(wǎng)格敏感性，于初始液位150 mm 工況下加密區(qū)域網(wǎng)格，分別采用了3 種不同的基礎尺寸，即0.006 m/0.003 m、0.005 m/0.002 5 m、0.004 m/0.00 2 m(采用較大尺寸標識)。圖4 為擋板擺角和最高液位對網(wǎng)格的敏感性。

圖4 最高液位與擋板擺角的網(wǎng)格敏感性

從擋板擺角來看，網(wǎng)格尺寸對模擬結果精度有一定影響。與0.006 m 網(wǎng)格相比，0.005 m 網(wǎng)格、0.004 m 網(wǎng)格結果有些許偏差；而0.004 m 網(wǎng)格結果顯示與0.005 m 網(wǎng)格結果相差無幾。從潰壩最高液位來看，網(wǎng)格尺寸越小，曲線越平滑。3 種網(wǎng)格尺寸情況下，網(wǎng)格越細越接近實驗值，最高液位對網(wǎng)格尺寸要求較高。綜合考慮精度和計算效率，模擬取加密區(qū)域0.005 m 網(wǎng)格尺寸計算較為合適。

4.2 計算結果驗證

擋板潰壩實驗中，潰壩水流泄出箱體，水流靜壓力作用在擋板上，擋板受到潰壩水流沖擊產生旋轉，在脫離沖擊后，擋板受重力作用開始做回復運動，在回復至一定角度受到潰壩水流泄出形成的潰壩水流沖擊，如此周期往復，如圖4(b)中擺角運動曲線。以下將以第一、二周期的方式來表述振蕩峰所在周期。

為了驗證模擬的有效性，需要將實驗結果與模擬結果進行比對。圖5 為在初始液位150 mm工況下，t=0.2 s 時，實驗和模擬在同一時間節(jié)點下的Y軸方向視角液位曲面形狀與擋板運動狀態(tài)的比對結果。圖6 為在初始液位150 mm 工況下，t=0.2 s 時，-X軸方向視角液位曲面形狀的實驗和模擬結果。圖5 和圖6 中能明顯觀察到潰壩水流泄出時鏡面水舌[11]兩側在前視窗上產生的液膜和完全潤濕區(qū)域。

圖5 Y 軸方向視角可視化結果示意

圖6 -X 軸方向視角可視化結果示意

圖7 為擋板角速度曲線，圖8 為初始液位150 mm 工況下，擋板受力、力矩圖。由圖7 和圖8 可以清晰地觀測到受力和力矩曲線有4 個峰，表明潰壩擋板受流體的影響具有周期性，且逐漸減弱。

圖7 擋板角速度

圖8 擋板所受力、力矩

圖9 為受力曲線峰值處擋板內表面壓力云圖。圖9 中：第一峰值處壓力大、受力面積大、壓力云圖梯度邊緣穩(wěn)定平滑；第二、第三峰值處壓力較小、受力面較小、梯度邊緣出現(xiàn)紊亂；第四峰值處所受力和受力面積均進一步減小。

圖9 各峰值時刻擋板內表面總壓力

4.3 最高液位變化分析

對擋板開啟運動過程中水箱最高液位變化進行分析，將液相體積份額為0.95 等值面的最高點作為最高液位，得到數(shù)值模擬最高液位數(shù)據(jù)。液面變化的實驗結果與模擬結果比對如圖10 所示。由圖10 可以看到總體趨勢模擬結果與實驗數(shù)據(jù)符合較好，說明VOF 方法適用于潰壩問題自由界面的模擬，本文模擬結果可信有效。在初始液位高度為100、150、200 mm 工況下，前期模擬值與實驗值符合較好；在約0.5 s 后，模擬值開始與實驗值產生些許偏差。

圖10 3 種工況下液面高度實驗值與模擬值

圖10 模擬結果呈階梯狀，這可能是由于潰壩水流經(jīng)擋板開口泄出形成的體積變化傳導至最高液位變化過程較緩，響應較慢，無法形成較平滑的變化曲線；另一方面，由于求解區(qū)域網(wǎng)格化的影響，使最高液位變化曲線呈階梯狀分布。結合圖4 網(wǎng)格敏感性分析結果，這說明VOF 方法對于時間響應較慢和自由界面捕捉情況下，對網(wǎng)格尺寸和時間步長要求較高，除非網(wǎng)格夠細或時間步長夠長，否則潰壩最高液面曲線無法平滑。

4.4 擋板開啟角度變化分析

圖11 對比了擋板開啟角度實驗結果與模擬結果，模擬結果反映了擋板運動角度變化趨勢，與實驗結果符合較好。

圖11 3 種工況下開啟角度實驗值與模擬值

擋板受潰壩水流形成的射流沖擊，其擺動呈現(xiàn)一定的規(guī)律性。如圖12 所示，隨初始液位的線性增加，擋板振動幅度呈遞增態(tài)勢，且遞增幅度逐漸減小。這是因為初始最高液位逐漸接近擋板旋轉軸，流體靜壓提供擋板的轉矩并非滿足線性遞增關系。振蕩周期越靠前，初始液位對振蕩幅值的影響越小，幅值遞增曲線越緩。

圖12 各周期振動幅值

隨初始液位的增加，周期時長呈遞減趨勢，如圖13 所示。值得注意的是，第二、第三周期時長較為接近，但第四周期時長卻遠大于第二、第三周期。這是由于第四周期處于擋板運動末期，擋板所受動壓較小，做功功率低，引起擋板運動狀態(tài)變化需要更長時間。

圖13 各周期時長

在初始液位遞增條件下，振蕩周期發(fā)生時間出現(xiàn)了明顯時延現(xiàn)象，如圖14 所示。振蕩周期越靠前，時延現(xiàn)象越不明顯，第一周期幾乎不會隨初始液位變化而時延。

圖14 各周期峰值時間點

與圖10 中潰壩最高液位曲線一樣，注意到初始液位100、150、200 mm 工況下，擋板擺角模擬值也均于約0.5 s 后開始偏離實驗值。結合圖11，偏差從擋板第一次回落受到射流沖擊時刻開始，根據(jù)圖9 壓力云圖峰值處出現(xiàn)壓力梯度邊緣不清晰現(xiàn)象，推測可能為VOF 模擬沖擊射流界面不夠精細導致DFBI 轉矩計算產生偏差，說明VOF 結合DFBI 模型對于潰壩水流形成的沖擊射流沖擊剛體運動還有進一步提升空間。

5 結論

1)本文通過將模擬與實驗結果對比，證明了VOF 方法結合DFBI 模型的CFD 數(shù)值模擬方法可以較精確地描述潰壩沖擊問題，為今后使用該數(shù)值模擬方法求解類似問題提供了有力的實驗依據(jù)。

2)潰壩擋板振蕩呈現(xiàn)周期性。振蕩幅度隨著初始液位的增加而增加，且隨初始液位增高，振蕩幅值增加程度逐漸降低，周期時長逐漸減小，時延現(xiàn)象逐漸減輕。

3)使用VOF 方法模擬潰壩水流自由表面，在時間響應較長和界面捕捉方面對于網(wǎng)格尺寸要求較高，適當加密網(wǎng)格會對結果大大改善。

4)VOF 方法和DFBI 模型結合對潰壩水流形成的沖擊射流沖擊剛體運動還有一定提升空間。