基于格子Boltzmann方法的潛堤上波浪破碎數(shù)值模擬

王藝之，韓新宇，董 勝

(中國海洋大學(xué)工程學(xué)院，山東 青島 266100)

引 言

近海海域的波浪，在行進中遇到潛堤時會發(fā)生劇烈破碎，隨后產(chǎn)生波浪增水，形成波生流以及低頻波運動。這些水動力過程會影響泥沙輸運、潛堤穩(wěn)定，甚至引起生物棲息地的變化，因此有必要對波浪在潛堤上破碎的情況進行全面的研究。

自19世紀50年代起，有學(xué)者針對潛堤上的波浪變形進行了觀察及試驗[1-5]。研究結(jié)果表明，陡坡與緩坡上的波浪變形非常相似，兩者的差異主要是波浪破碎的強度和位置。在陡坡上，水深沿程迅速變化，波浪劇烈破碎并產(chǎn)生一段狹窄的破波帶，水流可以看作是間歇性的強向岸射流后緊跟著緩慢的向?；亓鳌Ｅc此同時，針對潛堤上的波浪傳播問題，Symonds等[6]基于對自然礁上波浪的觀測以及對其能譜的分析，提出入射波的強度與淹沒水深有關(guān)；Hearn[7]研究了波浪驅(qū)動流在珊瑚礁上的流體動力學(xué)及其對礁體的沖刷；Gourlay和Colleter[8]利用明渠理論推導(dǎo)了礁頂流速與入射波條件、淹沒水深、堤頂物理特征的關(guān)系；Massel和Gourlay[9]提出，假設(shè)波和流的時間與空間尺度分離，與輻射應(yīng)力模型相結(jié)合的緩坡方程也適用于研究此類課題。

波浪在潛堤上破碎的數(shù)值模擬的關(guān)鍵，在于如何準確對復(fù)雜自由表面的變形進行建模。傳統(tǒng)的數(shù)值模擬方法是從宏觀得到微分方程，對該微分方程進行離散化計算并求解，只關(guān)心微分方程的形式與結(jié)果。而格子Boltzmann方法(LBM)并不是對方程直接求解，其基本思想是從介觀動力學(xué)出發(fā)，將流體離散成介觀粒子，按照符合物理規(guī)律的演化機制計算粒子系統(tǒng)，觀察流場的宏觀變化。該方法在宏觀上連續(xù)，微觀上離散，且邊界條件易處理，適合大規(guī)模的粘性流場計算[10]。

國內(nèi)最早由朱照宣[11]對LBM的前身格子氣自動機方法進行了總結(jié)，并提出了一些應(yīng)用方向；胡守信[12]于1989年在朱照宣[11]和U.Frisch[13]等的基礎(chǔ)上研究了直管道中繞平板平面流動的格子氣仿真，并于1997年與施衛(wèi)平共同研究了用LBM模擬無粘性淺水長波方程[14]是國內(nèi)較早對于LBM及其應(yīng)用的研究。隨后，LBM不斷發(fā)展，在三維風(fēng)場[15]、懸浮顆粒運動[16]、明渠水流模擬[17]等方面均得到了廣泛應(yīng)用。

本文探討LBM模型求解波浪在潛堤上破碎的適用性與準確性。數(shù)值建模時，采用推波板造波，并利用斜坡進行消波，構(gòu)造了二維的數(shù)值水槽。分別對高水位和低水位條件下的自由液面高度、波高的空間分布以及斜坡上波浪的增水進行了模擬，并將模擬結(jié)果與Wen等[18]的物理試驗數(shù)據(jù)進行了對比，此外還對不同工況下的波浪透射與反射系數(shù)進行了計算，所得結(jié)果對潛堤結(jié)構(gòu)設(shè)計有參考價值。

1 LBM數(shù)值模型

目前，LBM已經(jīng)逐漸成為解決流體流動方程的重要方法。該方法由元胞自動機、格子氣自動機模型結(jié)合Boltzmann方程發(fā)展而來。由于將格子類方法對速度、時間、空間的離散方式與Boltzmann方程的碰撞演化機制相結(jié)合，兼具二者的優(yōu)勢。

1.1 控制方程

采用BGK碰撞機制的Boltzmann方程表示如下：

(1)

式中，f=f(x,ξ,t)表示單個粒子的分布函數(shù)；ξ表示微觀速度；λ表示與碰撞有關(guān)的松弛時間；feq表示Boltzmann-Maxwellian平衡分布函數(shù)，該函數(shù)計算公式表示如下：

(2)

式中，D表示空間維數(shù)；ρ，u，T分別表示宏觀密度、速度、溫度。

將Boltzmann方程對時間進行離散，記沿特征線ξ的時間導(dǎo)數(shù)為:

(3)

把式(1)寫為關(guān)于時間的全微分形式，即:

(4)

對上式進行推導(dǎo)和泰勒展開，忽略高階無窮小項，記無量綱的松弛時間為τ=λ/Δt，可得對時間的離散結(jié)果如下：

(5)

在此基礎(chǔ)上，記{eα}為微觀速度ξ離散后的速度集合，則得對速度的離散結(jié)果為：

(6)

此即BGK逼近下的Boltzmann方程。

對于式(2)表示的分布函數(shù)，將其對速度進行泰勒展開至u2并截斷高階項，可得格子Boltzmann模型中的平衡分布函數(shù)為：

(7)

(8)

可以看出，選定離散速度{eα}后，只要再確定權(quán)系數(shù){ωα}就可確定具體的平衡分布函數(shù)。

對于n維空間中m個離散速度的情況錢躍竑提出DnQm系列模型，并隨LBM的發(fā)展不斷擴充，其中二維問題中最常用的就是D2Q9，也是本文模擬中所應(yīng)用的離散速度模型。該模型下的格子速度、權(quán)函數(shù)、當(dāng)?shù)芈曀偃缦拢?/p>

(9)

(10)

在該模型下通過多尺度展開可近似得到不可壓縮的N-S方程：

(11)

1.2 湍流模型

湍流模型采用Wall-Adapting Local Eddy-viscosity(WALE)模型，此模型對近壁面和遠壁面、層流和紊流都具有良好的特性[19]。該模型可以表示如下：

(12)

式中的系數(shù)定義如下：

(13)

式中，WALE常數(shù)Cw通常取0.2。

1.3 邊界條件

本文中應(yīng)用壁面函數(shù)如式(14)所示，同時考慮了曲率和壓力梯度的影響，能夠較準確的解決湍流邊界層的問題[19]：

(14)

式中，y+表示距離墻壁的法向距離；uτ表示表面摩擦速度；τw表示湍流壁面剪切應(yīng)力；U表示距離墻壁給定距離處的速度。具體定義如下：

(15)

2 數(shù)值模擬

數(shù)值模型設(shè)置如圖1所示。水槽長度由入射波長決定。水槽中放置潛堤模型，其陡坡坡度為1：1，坡后平臺長8.6 m，高1.8 m。在水槽起始端設(shè)置一推板，考慮減少推波板造成的二次反射，將推波板到潛堤的距離設(shè)為5.0L(L代表入射波波長)；在尾端設(shè)置消浪斜坡，坡度為1：10。水槽初始水深設(shè)置為d= 1.8+hr，其中hr代表潛堤的淹沒水深。

圖1 數(shù)值水槽示意圖

沿數(shù)值水槽共布置13個波高測點，它們的相對位置如圖1所示。G1—G3位于深水區(qū)，G4—G6位于陡坡區(qū)，G7—G13位于水平平臺區(qū)。坐標原點位于坡腳，傳感器的坐標如表1所示。

表1 波高儀坐標

考慮到計算精度以及計算效率，本文整體選擇4.0 cm的粒子間距，僅在自由表面以及潛堤附近使用粒子自適應(yīng)方法細化至2.0 cm，示意圖如圖2所示。模擬時長選取50 s。

圖2 自適應(yīng)粒子細化示意圖

3 基于數(shù)值結(jié)果的分析和討論

3.1 自由液面高度與試驗數(shù)據(jù)的對比

模擬自由液面高度時共設(shè)計兩個工況，波高均為0.20 m，周期均為1.75 s，淹沒水深分別為0.10 m(高水位)和0.02 m(低水位)。模擬結(jié)果如圖3和圖4所示。對比中用到試驗數(shù)據(jù)以及試驗圖片來自Wen的論文[18]。

圖3給出了兩個工況下試驗值和LBM計算值的對比結(jié)果。可以看出，計算值和試驗值吻合良好。說明該模型可以較為準確的模擬波浪在潛堤上的傳播過程。通過潛堤平臺上的測點結(jié)果，可以看到波浪在潛堤上傳播的非線性特征，波峰較陡，波谷較緩。由于潛堤的阻礙作用，平臺上的波高明顯減小。

圖3 試驗結(jié)果與LBM模擬結(jié)果對比

圖4是低水位下波浪與潛堤相互作用過程的自由表面的演變過程，左側(cè)為試驗圖，右側(cè)為數(shù)模圖。通過試驗和數(shù)值模擬可以看出，由于潛堤導(dǎo)致波面升高，隨后波浪在靠近斜坡的平臺處發(fā)生破碎。如圖4所示，當(dāng)t=0.0T時，波浪到達潛堤并開始破碎；當(dāng)t=0.25T時，崩破波在水平平臺傳播時，波浪表面開始出現(xiàn)液體飛濺；當(dāng)t=0.5T和t=0.75T時，已經(jīng)破碎的波浪繼續(xù)沿平臺傳播，與此同時，陡坡與平臺的交接處會發(fā)生輕微的回流。

圖4 瞬時波面對比(低水位)

3.2 波高的空間分布

模擬自由液面高度時，增加1.5 s、2.0 s兩種周期，共組合出6種工況進行數(shù)值模擬。波高空間分布的試驗與模擬對比情況如圖5所示，數(shù)值模擬沿程波高變化曲線由Akima插值得到，插值節(jié)點橫坐標與試驗結(jié)果相同，均為沿程測點的橫坐標，縱坐標為模擬結(jié)果在該測點位置的波高。從圖中可以看出，模擬結(jié)果與試驗基本吻合，波高隨著波浪的傳播方向整體呈逐漸衰減的趨勢。由于大部分的入射波能量從陡坡反射回深水區(qū)，因此在潛堤前可以觀察到局部駐波。此外，波浪在潛堤上的破碎導(dǎo)致波能消散，使得波高在潛堤上也迅速衰減。

圖5 波高在潛堤上的空間分布

3.3 反射系數(shù)

圖6表示潛堤在不同工況下的反射系數(shù)與入射波周期的關(guān)系。在潛堤堤腳前1倍入射波長與1.2倍入射波長的位置分別設(shè)置一個測點，根據(jù)這兩個測點處的波高時程變化可以計算得出每個工況的反射系數(shù)如表2所示。由圖6可以看出，當(dāng)?shù)躺纤顬?.10 m時，潛堤反射系數(shù)隨入射波周期增大而逐漸增大，且周期與反射系數(shù)近似成正比；當(dāng)?shù)躺纤顬?.02 m時，反射系數(shù)隨入射波周期的增加成先增加后減小的趨勢。兩種水位的反射系數(shù)都在入射波周期為1.5 s時最小，且入射波周期相同的情況下，低水位工況下潛堤的反射系數(shù)要比高水位工況的更大。綜上所述，潛堤在低水位工況下具有較好的消浪效果。由于本文涉及的工況下潛堤透射系數(shù)過小，可以忽略，故本文不對該潛堤的透射系數(shù)進行研究。

圖6 反射系數(shù)變化趨勢

表2 各工況反射系數(shù)

4 結(jié)論

基于LBM，本文對規(guī)則波與潛堤的相互作用進行了研究，建立二維波浪水槽對其進行數(shù)值模擬。通過模擬值和試驗值進行對比分析，結(jié)論如下：

(1)基于LBM方法可以較好的模擬規(guī)則波與潛堤相互作用過程中的波浪破碎和液體飛濺等非線性現(xiàn)象，傳播過程中的波高與試驗吻合良好。

(2)通過對沿程波高進行分析，波高整體呈衰減趨勢。在波浪傳播到靠近斜坡的平臺附近，由于波浪反射和破碎，波高迅速衰減。

(3)在相同波高的情況下，高水位反射系數(shù)隨周期的增大而增大，低水位反射系數(shù)隨周期的增大先增大后減小。