基于格子Boltzmann方法的斜坡堤越浪數(shù)值模擬研究

李薪豐，張慶河，張金鳳，劉光威, 2

(1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300350; 2.清華大學 水沙科學與水利水電工程國家重點實驗室，北京 100084)

斜坡式海堤是海防工程中常見的結(jié)構(gòu)物，海堤發(fā)生破壞可能給人們的財產(chǎn)與生命安全帶來重大損失。越浪往往是造成海堤破壞的重要因素，因此海堤越浪研究具有重要意義。早期的越浪研究主要集中在物理模型試驗，隨著波浪理論和計算機技術(shù)的發(fā)展，越來越多的數(shù)值方法，如基于N-S方程的OpenFOAM、Fluent、Flow-3D等開源或商業(yè)軟件[1-3]，以及基于拉格朗日方法的光滑粒子流體動力學(SPH)模型[4-5]，被應(yīng)用于越浪過程的研究。

近年來，格子Boltzmann方法(lattice Boltzmann method，簡稱LBM)作為新型的流體數(shù)值模擬方法，憑借其算法簡單，計算效率高，并行計算可擴展性好，能方便處理復(fù)雜邊界條件等優(yōu)勢，在涉及流體模擬的很多領(lǐng)域得到了廣泛應(yīng)用[6]，但目前應(yīng)用LBM對波浪問題進行模擬的工作還只見于少量文獻[7-10]。Liu等[11]指出，LBM方法在波浪研究中未能得到廣泛應(yīng)用的原因在于，已有文獻在利用LBM模擬波浪時存在較為嚴重的波浪衰減問題。Liu等[12]進一步分析了波浪衰減原因，并提出了壓力修正格式，解決了波浪衰減問題，為LBM應(yīng)用于波浪模擬奠定了基礎(chǔ)。

為進一步驗證修正格式LBM模擬波浪和結(jié)構(gòu)物相互作用的能力，將采用劉光威[13]建立的二維數(shù)值波浪水槽，對光滑斜坡堤上規(guī)則波與不規(guī)則波越浪進行模擬，并將計算結(jié)果與已有試驗數(shù)據(jù)以及其他數(shù)值模型的模擬結(jié)果進行比較。

1 數(shù)值模型

對于斜坡堤越浪問題，參考坐標系如圖1所示，坐標原點O固定在入口與底面相交處，x軸正方向為波浪的傳播方向，y軸正方向為結(jié)構(gòu)物的高度方向，SWL表示海邊界靜水位。

圖1 斜坡堤越浪問題的參考坐標系

1.1 LBM方法

在LBM模型中，流體抽象為大量的粒子，這些粒子在簡單的網(wǎng)格上進行碰撞和遷移，根據(jù)反映粒子分布狀態(tài)的分布函數(shù)和宏觀物理量之間的關(guān)系即可獲得密度和速度等宏觀流動信息。LBM只描述粒子的統(tǒng)計行為，而不關(guān)心粒子的個體行為，也可以看做是Boltzmann方程的特殊離散形式。

傳統(tǒng)的LBM在模擬波浪時，采用單相流模型，只對水相進行模擬，通過自由表面邊界條件來表示氣相對水相的影響，為考慮外力項(如重力)需要在LBM模型中添加作用力項，正是由于自由表面邊界條件與作用力項的耦合誤差，造成非物理流動引發(fā)了波能衰減與數(shù)值振蕩。為消除耦合誤差，Liu等[12]將重力引入壓力梯度項，采用如下不帶作用力項的多松弛(MRT)模型的控制方程：

(1)

(2)

式中：ωα為各個離散速度方向的權(quán)系數(shù)，在D2Q9模型中，ω0為4/9，ω1～ω4為1/9，ω5～ω8為1/36。其對應(yīng)的離散速度如下：

(3)

通過下式可以從分布函數(shù)中得到與流體運動有關(guān)的宏觀信息，用于計算碰撞步中的平衡態(tài)分布函數(shù)：

(4)

Liu等[12]采用修正壓強格式后，將重力引入壓力梯度并將其從LB方程右側(cè)的作用力項中移除，當模型中外力項只包含重力時，方程中就不存在作用力項，進而避免了由于自由表面邊界條件與作用力項耦合而帶來的誤差。

1.2 邊界條件

模型采用主動吸收式速度入口造波方法[13]。速度入口造波方法只需要給出計算邊界處水質(zhì)點運動速度以及水位。首先根據(jù)正確的波浪理論計算出造波邊界處的理論水位ηT以及水質(zhì)點的理論運動速度uT。然后根據(jù)理論水位ηT判斷造波邊界處各個格點的體積分數(shù)φ：

(5)

式中：yc為格點的垂向坐標，h為靜水位，r為緩啟動系數(shù)，其表達式如式(6)所示：

(6)

式中：Tramp是設(shè)置的緩啟動時間。根據(jù)理論水位、水質(zhì)點理論運動速度以及格點垂向坐標計算邊界格點的速度：

uc=(φ×r)uT(yu,t)

(7)

式中：yu為計算格點流速時的垂向坐標，當格點位于理論水位下方時，使用格點垂向坐標yc計算流速；當格點位于理論水位上方時，使用靜水位h計算流速。

模型采用出流邊界消波方法[13]，邊界處采用帶垂向分布的出流流速：

(8)

(9)

將出流邊界消波方法應(yīng)用于速度入口邊界條件，則可以實現(xiàn)吸收反射波的效果，此時ηR=ηT-ηM。主動吸收速度入口造波邊界上格點的水平流速為：

(10)

由于LB方法在計算過程中只求解密度分布函數(shù)，實際應(yīng)用中往往將邊界上的速度或壓強等宏觀量作為已知條件，因此采用非平衡態(tài)外推數(shù)值格式來實現(xiàn)邊界上由宏觀量重構(gòu)密度分布函數(shù)的過程[12-13]。

1.3 自由表面模型

界面追蹤模型采用Thuerey[14]提出的基于VOF方法的單相自由表面LB模型[12-13]。根據(jù)離散單元的體積分數(shù)φ追蹤自由表面，即按照體積分數(shù)來定義氣相、界面和液相格點：

(11)

式中：b為單元的孔隙率，體積分數(shù)φ定義如下：

(12)

其中，Al(x，t)為x格點所在單元體內(nèi)流體的面積，Ac為單元的面積，m(x，t)為單元體內(nèi)液相流體質(zhì)量。該模型從分布函數(shù)的物理意義出發(fā)，由格點x遷移至相鄰格點的分布函數(shù)可被視為從格點流出的質(zhì)量，而由相鄰格點遷移至格點x的分布函數(shù)被視為流入格點單元的質(zhì)量。則從t時刻到t+1時刻格點α方向的質(zhì)量變化Δmα為：

(13)

(14)

則t+1時刻的格點單元液相流體質(zhì)量為：

(15)

有了新時刻的質(zhì)量之后，通過式(12)計算所有界面格點新時刻的體積分數(shù)，然后根據(jù)式(11)更新界面格點新時刻的狀態(tài)。

2 光滑斜坡堤上越浪的LBM模擬

2.1 光滑斜坡堤上規(guī)則波的越浪模擬

Saville[15]對坡度為1∶3和1∶1.5的光滑斜坡堤上規(guī)則波越浪進行了大量的試驗研究。Soliman[16]和閆科諦[1]分別利用基于求解RANS方程和VOF方法的二維破波數(shù)值模型(BWNM)、三維OpenFOAM模型對越浪過程進行了詳細的數(shù)模研究，并與Saville[15]的試驗數(shù)據(jù)進行了比較。Shao等[4]和Ni等[5]分別利用基于SPH方法的二維ISPH模型、二維DualSPHysics模型也進行了類似的數(shù)值試驗。這里采用基于LBM建立的二維數(shù)值波浪水槽對規(guī)則波作用下光滑不透水斜坡堤越浪進行模擬，并將計算結(jié)果與Saville[15]試驗值以及其他模型的模擬值進行比較。采用無量綱平均越浪量Q作為比較參數(shù)，其定義如式(16)所示：

(16)

其中，H0是深水波高；g為重力加速度；q為平均越浪量，m3/(m·s)。q的計算方法為：模擬過程中實時輸出越過斜坡堤堤頂豎直面的流體體積通量φ，即為每個時間步的堤頂瞬時越浪量q(t),對q(t)進行積分運算，再除以積分段時間間隔，即可得到平均越浪量q。

試驗海堤剖面如圖2所示，海堤設(shè)置在坡度為1∶10的緩坡上，其中dt、ds、Rc分別表示海邊界靜水位(SWL)以下水深、結(jié)構(gòu)趾部SWL以下水深、SWL以上結(jié)構(gòu)頂面(干舷)高度。文中選取海堤坡度為1∶3的16組算例進行了模擬，16個算例的深水波高均為1.0 m。按照原型比尺進行模擬，詳細的試驗參數(shù)和計算結(jié)果如表1所示，其中Ht為造波邊界處的波高值，T為周期。

圖2 規(guī)則波作用下斜坡堤斷面

表1 規(guī)則波越浪試驗條件和結(jié)果

在LBM計算中，造波邊界距離緩坡堤腳約4倍波長。同物模試驗中的統(tǒng)計方法，計算越浪量時去除前幾個波浪越浪的結(jié)果，取第4至第5個波浪越浪過程Q的平均值作為結(jié)果。以算例6為例，由圖3的計算結(jié)果可以看出，統(tǒng)計時間段內(nèi)波浪越浪量基本已達到準穩(wěn)態(tài)。

圖4為試驗值與各模型模擬值的比較，其中LBM的模擬結(jié)果整體上與試驗數(shù)據(jù)吻合較好，在合理的精度范圍內(nèi)，表明對越浪的預(yù)測沒有系統(tǒng)性偏差。最差的結(jié)果出現(xiàn)在算例8和算例11中，這兩組的越浪量試驗值在所有算例中最小，被嚴重高估。除了Soliman[16]的二維BWNM模型，其他模型在算例8和算例11中也均出現(xiàn)了越浪量被嚴重高估的現(xiàn)象，且均為偏差最大的組(如圖4所示)。Soliman[16]的模型中越浪量被嚴重高估的兩組為算例11和算例7，其中算例7同樣為越浪量試驗值較小的情況。

對于其他算例LBM計算的平均偏差為24.50%，同樣去除算例8和算例11的情況下Shao等[4]和Ni等[5]的SPH模型的平均偏差分別為18.88%和17.20%。閆科諦[1]采用OpenFOAM模型僅計算了算例1～11，故剩余9組的平均偏差為24.38%。Soliman[16]去掉誤差最大的兩組后平均偏差為20.24%。由以上數(shù)據(jù)可見，SPH模型的越浪量模擬誤差最小，各數(shù)值模型的預(yù)測水平基本相當，LBM的模擬結(jié)果略差于其他數(shù)值模型。

各模型在越浪量較小，越浪過程中破碎比較劇烈的情況下模擬偏差均較大，可能與實驗室在越浪較小時測量誤差較大有關(guān)[4-5]。LBM相較于其他模型誤差偏大的原因還在于二維條件下LES紊流模型不完全適用，未來可以通過引入RNG k-ω等模型來改善模擬結(jié)果。另外，與采用拉格朗日法的SPH模型相比，LBM、BWNM、OpenFOAM模型均采用歐拉網(wǎng)格法，使用的基于VOF方法的自由表面模型僅為一階精度，精度較低，在處理自由表面問題方面的性能比不上SPH模型，因此模擬誤差較大。未來可進一步在LBM中開發(fā)高精度的自由表面模型來彌補此不足。

圖3 算例6的各個周期的平均越浪量

圖4 規(guī)則波作用下不同模型越浪量計算值與試驗值比較

2.2 光滑斜坡堤上不規(guī)則波的越浪模擬

Van der Meer和Janssen[17]針對不規(guī)則波作用下海堤越浪進行了大量試驗，并進一步建立了越浪公式。Soliman[16]和Shao等[4]分別利用Van der Meer和Janssen[17]試驗數(shù)據(jù)對其二維BWNW模型、ISPH模型進行了進一步驗證。這里也將使用二維LBM模型對不規(guī)則波作用下坡度為1∶1、1∶2和1∶4的光滑不透水斜坡堤越浪進行模擬，并將計算結(jié)果與Van der Meer和Janssen[17]試驗值及其他模型模擬值進行比較。

試驗斷面如圖5所示，各算例靜水深dt均為8.0 m，其中Rc為干舷高度，造波邊界距離斜坡堤堤腳約3倍有效波長。采用Goda[18]修正的Jonswap譜生成不規(guī)則波，波譜增強因子γ=3.3。在計算中，頻譜由50個分量頻率表示，有效頻率范圍為0.5fp至3.5fp，fp為譜峰頻率。各算例詳細參數(shù)如表2所示，其中Hs為有效波高，Tp為譜峰周期，Tm為平均波周期。

圖5 不規(guī)則波作用下斜坡堤斷面

表2 不規(guī)則波越浪試驗條件和結(jié)果

算例4的累計越浪量歷時曲線如圖6所示，從圖中可以看出不規(guī)則波越浪表現(xiàn)出明顯的隨機性，各個周期的越浪量明顯不同，因此采用平均越浪量q作為評價指標研究不規(guī)則波越浪。根據(jù)Soliman[16]的研究，模型運行30 s左右波浪越浪量趨于穩(wěn)定。Shao等[4]在SPH的計算中發(fā)現(xiàn)越浪需要50 s的時間才能變得穩(wěn)定，故采用50～180 s內(nèi)的模擬結(jié)果計算平均越浪量。在LBM計算中，發(fā)現(xiàn)需要15Tp的時間越浪過程才能穩(wěn)定，即66.45～91.20 s的時間。圖7為從15Tp開始統(tǒng)計的算例4的平均越浪量歷時曲線，可以看出255Tp以后平均越浪量開始趨于穩(wěn)定，因此文中取15Tp～315Tp時間段內(nèi)的結(jié)果來計算平均越浪量。

圖6 算例4的累計越浪量歷時曲線

圖7 算例4的平均越浪量歷時曲線

圖8顯示了LBM的模擬結(jié)果與Van der Meer和Janssen[17]實驗室數(shù)據(jù)，以及Soliman[16]利用BWNM模型模擬結(jié)果和Shao等[4]利用ISPH模型模擬結(jié)果的對比。結(jié)果表明，各個模型均能較好地預(yù)測平均越浪量。LBM的預(yù)測平均偏差為19.31%，相比之下，BWNM模型的預(yù)測平均偏差為27.77%，LBM模型的模擬結(jié)果略優(yōu)于Soliman[16]的模型。SPH模型僅計算了坡度為1∶1和1∶4的算例，故算例1～4與算例10～11的平均偏差為15.44%，預(yù)測的準確性方面在三個模型中表現(xiàn)最好。

根據(jù)Soliman[16]的描述，越陡的海堤越容易產(chǎn)生非破波越浪，而越平緩的海堤越容易產(chǎn)生破波越浪。這與圖9和圖10所示的坡度分別為1∶1、1∶4的算例4、算例10的越浪過程波面圖相一致。LBM對于坡度為1∶1、1∶2、1∶4算例的平均偏差分別為11.18%、12.52%、52.53%，由此可以看出LBM對于破碎情況越劇烈的越浪，模擬誤差越大。因此，可以進一步提高LBM自由表面模型的精度以及開發(fā)合理的二維紊流模型來改善破碎波模擬結(jié)果，進一步改善LBM模擬精度。

圖8 不規(guī)則波作用下不同模型越浪量計算值與試驗值比較

圖9 算例4的越浪過程波面圖

圖10 算例10的越浪過程波面圖

3 結(jié) 語

采用劉光威[13]基于LBM方法建立的二維數(shù)值波浪水槽，通過主動吸收式速度入口造波，出流邊界消波的方法對光滑斜坡堤在規(guī)則波與不規(guī)則波作用下的越浪進行了模擬，模擬結(jié)果與試驗值吻合較好，與N-S方程模型和SPH模型模擬精度基本相當，證明了此模型具有模擬斜坡堤越浪的能力。若要對小越浪量等波浪破碎劇烈的情況實現(xiàn)更為準確的模擬，LBM模型還需要進行不斷發(fā)展與完善，進一步提高自由表面模型的精度以及開發(fā)更為適用的二維紊流模型。

就計算時間來看，LBM具有較高的并行計算效率，與基于不可壓縮N-S方程的傳統(tǒng)CFD模型和SPH相比有更好的并行可擴展性[13]。以文中的不規(guī)則波算例為例，這些算例的模擬均在“天河1-A”超算平臺上進行，使用未加密的均勻正方形網(wǎng)格，如算例4模擬時間共1 915.2 s，使用48個CPU核，模擬網(wǎng)格數(shù)為91 540個，計算用時為30 min。LBM數(shù)值波浪水槽的建立有望為研究原型尺度的三維波浪和結(jié)構(gòu)物相互作用提供更為迅捷的模擬手段，值得進一步推廣。