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

(1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350; 2.清華大學(xué) 水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

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

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

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

1 數(shù)值模型

對(duì)于斜坡堤越浪問(wèn)題，參考坐標(biāo)系如圖1所示，坐標(biāo)原點(diǎn)O固定在入口與底面相交處，x軸正方向?yàn)椴ɡ说膫鞑シ较颍瑈軸正方向?yàn)榻Y(jié)構(gòu)物的高度方向，SWL表示海邊界靜水位。

圖1 斜坡堤越浪問(wèn)題的參考坐標(biāo)系

1.1 LBM方法

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

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

(1)

(2)

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

(3)

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

(4)

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

1.2 邊界條件

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

(5)

式中：yc為格點(diǎn)的垂向坐標(biāo)，h為靜水位，r為緩啟動(dòng)系數(shù)，其表達(dá)式如式(6)所示：

(6)

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

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

(7)

式中：yu為計(jì)算格點(diǎn)流速時(shí)的垂向坐標(biāo)，當(dāng)格點(diǎn)位于理論水位下方時(shí)，使用格點(diǎn)垂向坐標(biāo)yc計(jì)算流速；當(dāng)格點(diǎn)位于理論水位上方時(shí)，使用靜水位h計(jì)算流速。

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

(8)

(9)

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

(10)

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

1.3 自由表面模型

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

(11)

式中：b為單元的孔隙率，體積分?jǐn)?shù)φ定義如下：

(12)

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

(13)

(14)

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

(15)

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

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

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

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

(16)

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

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

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

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

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

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

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

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

圖3 算例6的各個(gè)周期的平均越浪量

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

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

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

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

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

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

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

圖6 算例4的累計(jì)越浪量歷時(shí)曲線

圖7 算例4的平均越浪量歷時(shí)曲線

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

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

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

圖9 算例4的越浪過(guò)程波面圖

圖10 算例10的越浪過(guò)程波面圖

3 結(jié) 語(yǔ)

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

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